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PREFACE 


This  work  was  performed  under  U.S.  Navy  Contract  N00014-83- 
C-0650.  PSI  employees  who  contributed  substantially  to  the 
content  of  this  report  include  Mr.  Wilburt  Geddes  who  assisted 
greatly  in  the  identification  of  geophysical  data  sources, 
Messrs.  Robert  Hagg,  David  Delaney,  and  Michael  Switney  who 
performed  the  manual  data  entry,  Mr.  Douglas  Waugh  who  generated 
the  software  for  digital  data  entry  and  the  Monte  Carlo 
Biot/Stoll  runs,  and  Mr.  Daniel  B.  Asher  who  generated,  operated 
or  modified  all  the  remaining  software.  The  holdings  of  the 
National  Geophysical  Data  Center,  Boulder,  Colorado  were  the 
primary  source  for  the  data  used  in  this  research.  The  REFLEC 
model  was  developed  by  the  Naval  Ocean  Research  and  Development 
Activity  under  the  Bottom  Interaction  Program  and  made  available 
for  this  effort.  The  Bottom  Interaction  Program  also  sponsored 
the  development  at  PSI  of  the  PHYSED  and  SETUP  software  version 
of  the  Biot/Stoll  model  which  provided  software  modules  used  in 
the  present  study.  In  addition,  the  Bottom  Interaction  Program 
sponsored  a  PSI  study  of  depth  dependence  of  geophysical 
properties  that  led  to  the  selection  of  the  procedure  used  in  the 
present  work. 
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.--Bottom  loss  (BL)  can  be  computed  from  profiles  of  the  geoacoustic  parameters  density, 
compressional  speed  and  compressional  attenuation.  The  geoacoustic  properties  can  be 
derived  from  physical  properties  of  the  sediments,  which  are  more  readily  available, 
using  a  physical  sediment  model  (PHYSED)  based  on  the  Biot  theory  of  acoustic 
propagation  in  porous  media.  We  predict  BL  for  sediments  on  the  continental  terrace 
and  evaluate  those  predictions  in  terms  of  their  precision.  We  use  readily  available 
physical  properties  assembled  into  the  data  base  PHYPR0SE,  with  the  PHYSED  model  to 
produce  geoacoustic  profiles  for  the  calculation  of  BL  by  the  computer  program  REFLEC. 
We  divide  continental  terrace  sediments  into  provinces  based  on  water  depth  ranges  and 
grain  size  classes.  We  find  that  useful  separations  of  BL  values  by  province  do 
occur.  Some  provinces  show  significantly  different  mean  BL  values  and  some  provinces 
show  significantly  reduced  variances.  This  result  has  particular  importance  in  the 
shallow  water  environment,  since  sediment  physical  properties  are  often  the  only 
readily  obtainable  information  describing  the  seafloor.  We  strongly  recommend  the 
application  of  this  physical  sediment  approach  in  shallow  water  to  further  the  devel¬ 
opment  of  seafloor  interaction  models  and  data  bases  suitable  for  operational  use. 

Anticipated  Benefits/Potential  Commercial  Applications  of  the  Research  or  Development 

The  PHYSED  model  and  the  PHYPR0SE  data  base  have  been  used  in  research  supporting 
sediment  province  determinations  for  mine  warfare  and  ASW  applications.  Their 
continued  use  in  such  applications  is  indicated  through  funding  by  the  Acoustic  Bottom 
Interaction  Program  at  N0RDA,  and  the  AEAS  Program  at  0NR. 


1.0  INTRODUCTION 


1 . 1  Background 

Navy  acoustic  ASW  system  performance  can  be  affected  by  the 
properties  of  the  seafloor.  The  seafloor  can  act  as  part  of  the 
acoustic  waveguide  that  connects  source  and  receiver,  as  a 
barrier  blocking  such  transmissions,  or  as  a  scattering  element 
contributing  to  the  variance  in  acoustic  system  performance  on 
short  scales.  The  weapons  and  sonar  systems  affected  by  bottom 
interaction  include  systems  designed  to  exploit  bottom-bounce 
paths  (SQS-26/53  surface  ship  sonars,  BQS-13  and  BQQ-5  Sphere 
submarine  sonars,  VLAD  air  deployed  sonobuoys),  all  systems 
operating  in  a  "close  and  localize"  mode,  emerging  broadband 
localization  systems  (WAA,  RAPLOC,  LSI  concepts,  ACSAS,  TARP) , 
active  adjunct,  shallow  water  torpedo  sonars,  mine  hunting 
sonars,  and  bottom  or  near-bottom  mounted  surveillance  systems. 
The  frequency  range  of  interest  for  performance  estimates  of  such 
systems  spans  from  less  than  100  Hz  to  greater  than  106  Hz.  For 
a  given  system  and  geometry,  whether  the  seafloor  acts  as  barrier 
or  waveguide  depends  upon  its  acoustic  properties.  This  report 
is  concerned  with  the  .problem  of  predicting  the  acoustic 
properties  of  the  seafloor  in  terms  useful  for  assessing  the 
performance  of  ASW  systems. 

The  measure  of  seafloor  interaction  employed  in  the  sonar 
equations  is  bottom  loss  (BL).  BL  is  expressed  in  decibels  (dB) 
relative  to  the  energy  of  the  incident  acoustic  wave  and  is  a 
function  of  frequency  and  grazing  angle.  BL  is  controlled  by 
compressional  wave  speed  and  attenuation,  shear  wave  speed  and 
attenuation,  and  sediment  density.  These  we  term  the  geoacoustic 
properties.  Methods  that  can  predict  bottom  loss  or  geoacoustic 
properties  for  a  given  location  are  extremely  valuable. 
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The  Navy  recently  developed  the  Bottom  Loss  Upgrade  (BLUG)-'- 
model  and  data  base  to  provide  bottom  loss  for  deep  ocean 
locations.  BLUG  defines,  for  a  location,  smoothed  vertical 
profiles  of  the  geoacoustic  properties  that  produce  BL  curves 
essentially  equivalent  to  measurements  of  BL  at  that  location. 
Locations  devoid  of  BL  measurements  are  associated  with  nearby 
geoacoustic  profiles  derived  from  BL  if  sediment  type  and 
thickness  are  similar.  BLUG  defines  a  single  curve  for  shallow 
water,  continental  shelf  and  slope  sediments.  Given  the  extreme 
variability  of  shallow  water  sediments  compared  to  deep  ocean 
sediments  and  the  strategic  importance  of  shallow  water  straits, 
coastal  boundaries,  and  choke  points,  a  single  shallow  water 
geoacoustic  profile  is  inappropriate  for  many  applications. 

Because  direct  measurements  of  the  geoacoustic  properties 
are  costly  and  time  consuming,  other  approaches  have  been 
developed.  One  approach,  developed  and  pursued  by  Hamilton2 3 4,  is 
to  establish  empirical  relationships  between  geoacoustic 
properties  of  a  sediment  and  physical  properties  which  are  more 
numerous  or  less  costly  to  accumulate.  Another  approach  is  to 
relate  the  geoacoustic  properties  to  the  physical  properties 
based  on  physical  principles  and  a  single  comprehensive  theory  of 
porous  media  developed  by  Biot2»^»^  and  applied  to  saturated 

^Spofford,  C.W.  ,  R.R.  Greene,  and  J.B.  Hersey,  1983,  "The 
Estimation  of  Geo-Acoustic  Ocean  Sediment  Parameters  from 
Measured  Bottom-loss  Data",  SAI -83-879-WA ,  Science  Applications 
International  Corp.,  McLean,  VA. 

2Hamilton,  E.L.,  1980,  "Geoacoustic  Modeling  of  the  Sea  Floor", 

J.  Acoust.  Soc .  Am..  68,  pp  1313-1340. 

3 

Biot,  M.A.,  1956  "Theory  of  Elastic  Wave  Propagation  in  a  Fluid- 
Saturated  Porous  Solid",  I.  Low  Frequency  Range,  J.  Acoust.  Soc. 
Am. .  28,  pp  168-178. 

4 

Biot,  M.A.,  1956,  "Theory  of  Elastic  Wave  Propagation  in  a 
Fluid-Saturated  Porous  Solid”,  II.  Higher  Frequency  Range,  J . 
Acoust .  Soc .  Am . .  28,  pp  179-191. 

CL 

Biot,  M.A.,  1962,  "Generalized  Theory  of  Acoustic  Propagation  in 
Porous  Dissipative  Media",  J .  Acoust .  Soc .  Am . ,  34,  pp  1254-1264 


marine  sediments  by  Stoll.  ’’  The  advantages  of  the  Biot/Stoll 
approach  over  the  Hamilton  and  BLUG  approaches  have  been 
described  elsewhere.  ’  They  include  allowing  a  more  realistic 
dependence  of  attenuation  on  frequency  for  a  broader  range  of 
sediments,  providing  internally  consistent  geoacoustic 
descriptions,  and  requiring  fewer  empirical  relationships  to 
generate  a  full  geoacoustic  description  (thereby  accelerating  the 
description  process).  The  accuracy  of  the  Biot/Stoll  approach 
(as  implemented  and  extended  in  the  computer  program  PHYSED)  in 
marine  sediments  has  been  demonstrated  to  be  as  good  as  the 
measurements  to  which  it  was  being  compared^’ when  the 
inputs  were  properly  defined. 

To  be  useful  as  a  predictor  of  BL  at  a  shallow  water  loca¬ 
tion,  the  physical  properties  upon  which  PHYSED  depends  must  be 
in  sufficient  supply  to  great  enough  precision  that  uncertainties 
in  the  resulting  BL  are  acceptably  small.  Improvement  over 
existing  capability  is  attained  if  the  precision  is  on  the  order 
of  that  for  BL  measurements,  or  better,  and  physical  property 
availability  is  greater  than  BL  and  geoacoustic  measurements  com¬ 
bined.  It  is  the  purpose  of  the  work  presented  here  to  estimate 
the  availability  and  the  precision  of  BL  predictions  based  upon 


Stoll,  R.D. ,  1974-  "Acoustic  Waves  in  Saturated  Sediments",  in 
L.  Hampton  (Ed.),  Physics  of  Sound  in  Marine  Sediments. 


Stoll,  R.D.,  1974,  "Acoustic  Waves  in  Ocean  Sediments", 
Geophysics .  42,  pp  715-715. 

Stoll,  R.D.,  1980,  "Theoretical  Aspects  of  Sound  Transmission  in 
Sediments",  J.  Acoust.  Soc .  Am..  68,  pp  1341-1350. 

Brunson,  B.A.,  and  E.J.  Molinelli,  1982,  "A  Physical  Sediment 
Model  for  the  Prediction  of  Seafloor  Geoacoustic  Properties", 

PSI,  TR-216227  for  ONR ,  Planning  Systems  Inc.,  McLean,  VA,  22101 

^Holland,  C.W. ,  and  B.A.  Brunson,  1985,  "The  Biot/Stoll  Model: 

An  Experimental  Assessment”,  PSI  TR-185331  for  NORDA  Code  113, 
Planning  Systems,  Inc.,  McLean,  VA.,  22102 

^Beebe,  J.H.,  198'>,  "An  Experimental  Investigation  of  Ocean 
Sediment  Effect.  Upon  Long-Range  Transmission  Loss  in  Shallow 
Water",  Technical  Memorandum  TM  80-247,  Pennsylvania  State  Univ. 


the  PHYSED  model.  By  precision  we  mean  the  spread,  or  variance, 
or  uncertainty,  in  BL  values  associated  with  a  given  prediction. 


The  PHYSED  model  is  a  formulation  that  relates  the  sound 
speed  and  attenuation  of  acoustic  waves  in  the  sediment  to  the 
physical  properties  of  sediment  constituents  --  the  sediment 


grains,  the  pore  fluid,  and  the  structure  of  the  grains  within 
the  sediment  (the  dry  "frame").  Table  1-1  lists  the  sediment 
constituent  properties  that  are  required  as  input  to  the  PHYSED 
model.  These  properties,  together  with  their  units  and  typical 
values,  are  discussed  in  some  detail  by  Brunson  and  Molinelli.® 


Table  1-1.  PHYSED  Model  Physical  Inputs 


SYMBOL 

PHYSICAL  PROPERTY 

UNITS 

pr 

Grain  Properties 

Density  of  sediment  grains 

g  cm-3 

K 

r 

Bulk  Modulus  of  sediments  grains 

-2 

dyne  cm 

Pf 

Pore  Fluid  Properties 

Density  of  pore  fluid 

g  cm-3 

Kf 

Bulk  modulus  of  pore  fluid 

dyne  cm-2 

n 

Viscosity  of  pore  fluid 

g  cm-^  s — 1 

6 

Frame  Properties 

Porosity 

k 

Permeability 

_2 

cm 

a 

Pore  size  parameter 

cm 

a 

Structure  factor 

-- 

pb 

Shear  modulus  of  frame  (real  part) 

.  -2 
dyne  cm 

pb* 

Shear  modulus  of  frame 

_  o 

dyne  cm 

Kb 

(imaginary  part) 

Bulk  modulus  of  frame  (real  part) 

,  -2 
dyne  cm 

V 

Bulk  modulus  of  frame 

_  p 

dyne  cm 

(imaginary  part) 

vvy. 
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The  output  of  the  PHYSED  model  is  the  complete  frequency- 
dependent  description  of  compressional  wave  speeds  and 
attenuations  and  shear  wave  speeds  and  attenuation  valid  in  the 
frequency  range  from  10  to  10®  Hertz.  The  output  is  listed  in 
Table  1-2.  The  type  II  compressional  wave  is  a  diffusion  type 
wave  that  is  only  noticeably  excited  in  materials  in  which  the 
interstitial  fluid  has  a  small  bulk  modulus  compared  to  the 
frame.  It  is  not  important  in  most  marine  sediment  applications. 
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Table  1-2.  PHYSED  Geoacoustic  Output  Parameters 


Symbol 


Geoacoustic  Property 


Units 


Compressional  wave  speed  -  Type  I  m  s” 

Compressional  wave  attenuation  -  Type  I  dB  m-1 

Compressional  wave  speed  -  Type  II  ms-1 

Compressional  wave  attenuation  -  Type  II  dB  m-1 

_  .  _i 

Shear  wave  speed  m  s 

Shear  wave  attenuation  dB  m-1 


To  use  the  PHYSED  model  to  derive  the  geoacoustic  properties 
at  a  location  it  is  necessary  to  provide  all  thirteen  inDuts 
listed  in  Table  1-1.  However,  for  a  given  location,  values  for 
all  thirteen  physical  properties  are  rarely,  if  ever,  available. 
This  makes  it  necessary  to  obtain  several  of  the  thirteen  inputs 
by  other  means.  Some  of  them  may  be  derived  from  other  available 
properties  using  relationships  based  on  physical  principles. 
Others  may  be  derived  using  emp i r i ca 1  relationships  to  available 
properties.  Finally,  some  may  have  to  be  assigned  values  within 
some  global  or  local  range  known  to  apply  for  them.  The  specific 
options  so  far  assembled  for  obtaining  the  thirteen  inputs  have 


been  described  in  reference 


The  combined  use  of 


empirical  relationships  as  introduced  by  Stoll  and  others  with 
the  theory  of  Biot  has  led  to  our  use  of  the  term  "hybrid  model" 
for  the  suite  of  formulas  assembled  as  the  model  PHYSED. 12 » 13 » 14 
Figure  1-1  diagrams  the  relationship  between  the  thirteen  inputs 
(along  the  bottom  of  the  figure)  and  these  other  properties, 
which  we  term  the  "genesis"  parameters,  along  the  left  margin  of 
the  figure.  Table  1-3  lists  the  genesis  parameters,  their  sym¬ 
bols  and  units.  Computer  software  that  enables  a  user  to  run  the 
model  on  either  the  input  or  genesis  parameters  was  developed  as 
the  program  PHYSED.  This  software  had  to  be  modified  for  the 
purposes  of  the  present  effort  as  described  in  our  approach. 

The  central  issue  in  assessing  the  operational  utility  of 
PHYSED  modeling  is  whether  the  input  or  genesis  parameters  can  be 
provided  with  sufficient  horizontal  and  vertical  resolution,  with 
sufficient  global  coverage,  with  enough  accuracy,  with  the  needed 
precision  to  obtain  useful  geoacoustic  descriptions  of  the 
seafloor.  Usefulness  implies  that  the  model  inputs  are  more 
easily  obtained  than  collecting  and  compiling  measurements  of  BL 
and  geoacoustic  profiles.  The  first  part  of  the  assessment 
requires  assembling  a  data  base  of  available  physical  properties 
and  manipulating  them  so  as  to  provide  complete  inputs  for  the 
model.  If  enough  information  is  available  to  drive  the  PHYSED 
model  then  the  assessment  requires  that  the  observed  variability 
of  the  inputs  be  carried  through  to  the  variability  of 
geoacoustic  outputs  and  thence  to  the  variability  in  predicted 
BL.  PHYSED  modeling  will  be  considered  operationally  useful  if 

■'"^Ogushwitz,  P.R.,  1983,  "Appl icabl i ty  of  the  Biot  Theory,  I.  Low 
Porosity  Materials",  J.  Acoust.  Soc .  Am,. 

l^Molinelli,  E.J.,  and  B.A.  Brunson,  1983,  "PHYSED:  Physical 
Sediment  Model  Software  --  Technical  Manual",  PSI  TR-185246 
for  NORDA  Code  113,  Planning  Systems,  Inc.,  McLean,  VA,  22102. 

■^Watters,  P.D.,  1983,  "PHYSED:  Physical  Sediment  Model  Software 
—  User's  Manual",  TR-185245  for  NORDA  Code  113,  Planning 
Systems  Inc.,  McLean,  VA,  22102. 
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Table  1-3.  Genesis  Physical  Parameters 
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Symbol 

Physical  Property 

Units 

Z 

Depth 

m 

- 

Size  class  (sand,  silt,  clay) 

- 

- 

Mineral  content  (calcareous,  etc.) 

- 

dmg 

Mean  grain  diameter 

cm 

4> 

-log2  of  grain  diameter  in  mm 

phi  units 

<J<j> 

Standard  deviation  of  $  value 

phi  units 

So 

Specific  surface  (area/volume) 

cm-1 

Pr 

Shear  modulus  of  grains 

dynes  cm-2 

% 

Poisson's  ratio  of  sediment  frame 

- 

V 

s 

Shear  wave  speed  in  frame 

-1 

cm  sec 

VE 

Longitudinal  wave  speed  in  frame 

-1 

cm  sec 

As 

Frame  shear  wave  log  dec 

- 

ae 

Frame  longitudinal  wave  log  dec 

- 

- 

Inclusion  geometry 

- 

'  AP 

Frame  compressional  wave  log  dec 
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the  variability  in  predicted  BL  is  no  worse  than  what  can  now  be 
associated  with  a  "province",  whether  by  theory  or  measurements. 
We  will  define  a  "province"  later  in  the  report. 

Table  1-4  illustrates  the  stages  of  development  of  the 
PHYSED  model.  The  model  is  considered  at  stage  IV  by  virtue  of 
the  results  of  Holland  and  Brunson.10  The  objective  of  the 
present  study  is  to  develop  and  evaluate  the  PHYSED  model  for  use 
in  the  predictive  mode.  This  report  therefore  records  the 
progress  of  the  model  through  stage  V  and  and  into  stage  VI.  The 
physics  of  the  formulation  having  been  tested  and  the  utility  of 
the  hybrid  model  in  field  applications  having  been  shown,  we  now 
address  the  question  of  how  well  BL  can  be  predicted  for  a  given 
location  in  shallow  water. 

1.2  Objective 

Here  we  state  our  objective  again  for  clarity  and  emphasis. 
The  objective  of  these  efforts  is  to  estimate  how  well  we  can 
predict  BL  at  a  given  shallow  water  site.  We  will  use  four 
criteria  to  assess  the  quality  of  our  predictions. 

1)  The  precision  of  the  BL  value  assigned  to  a  class,  or 
province,  of  sediments  --  expressed  in  terms  of 
estimates  of  BL  variance  in  dB. 

2)  The  location  of  the  BL  distribution  for  a  class  -- 
expressed  in  terms  of  the  estimated  mean  BL  in  dB 

3)  The  accuracy  of  the  geoacoustic  values  associated  with 
a  province  --  expressed,  if  possible,  in  terms  of 
estimates  of  probability  that  a  given  measurement  of  BL 
or  a  geoacoustic  property  is  not  a  member  of  the 
distribution  assigned  to  the  province 


Table  1-4.  Stages  in  the  Development  of  a  Geoacoustic 
Description  Using  the  Biot/Stoll  Approach 


Stage 


Name 


Theoretical  Model 
Development 


Theoretical  Model 
Validation 


Hybrid  Model 
Development 


Hybrid  Model 
Validation 


Predictive  Model 
Development 


Predictive  Model 
Evaluation 


Predictive  Model 
Validation 


Research  Topics 


Principles  and  Mechanisms 


S-SS? 


■ 


Measurement/model  comparisons 
simple,  laboratory  sediments 

Measurement/model  comparisons 
in  complex,  laboratory 
sediments 


Published  measurements  and 
relationships  and  literature 
searches 


Empirical  relationship 
adaptation 


Measurement/model  comparisons 
for  in  situ  sediments 


Assemble  a  data  base  of  hybrid 
model  inputs 


Determine  precision  of  inputs 
Determine  precision  of  outputs 


Couple  geoacoustic  description 
to  an  acoustic  model 


Determine  effect  of  geo¬ 
acoustic  variations  on 
acoustic  model  outputs 

Compare  acoustic  model 
outputs  to  acoustic 
measurements 


Extend  data  base  of  hybrid 
model  inputs 

Design  and  perform  EVA 
experiment(s) 

Validation  experiment/model 
comparisons 

Operational  evaluation 


4)  the  availability  of  information  necessary  to  assign  a 
geographical  location  to  a  province  --  expressed  in 

terms  of  estimates  of  percent  of  available  information. 

We  consider  that  for  the  Biot/Stoll  model  to  represent 
potential  for  real  improvement  in  capability  the  information 
needed  to  assign  a  location  to  a  province  should  be  at  least 
twice  as  available  as  direct  measurements  of  geoacoustic  profiles 
and  BL  combined.  In  addition  at  least  one  of  the  following  must 
occur . 

1)  The  precision  found  for  a  province  must  be  greater 

(i.e.  have  significantly  lower  variance)  than  that  for 
shallow  water  sediments  as  a  whole. 

2)  The  mean  BL  for  a  province  must  be  significantly 

different  than  the  mean  BL  for  shallow  water  sediments 
overall . 

3)  The  probability  that  a  measurement  is  not  from  the 

population  ascribed  to  its  province  should  be  less  than 
one  percent. 

1.3  Approach 

Our  basic  approach  is  to  use  real  data  to  characterize  the 
variance  (uncertainty)  of  the  physical  properties  that  drive  the 
PHYSED  model,  to  propagate  that  variance  through  to  geoacoustic 
profiles  using  PHYSED,  to  propagate  that  variance  further  to  BL 
with  the  REFLEC  model,  and  to  describe  the  resulting  precision. 
The  details  of  the  approach  require  many  steps,  some 
approximations,  and  a  few  assumptions.  These  are  described  in 
the  remainder  of  this  section. 

The  first  step  is  to  assemble  data  into  a  computer  data 
base.  For  this  effort  data  sources  have  to  be  identified,  a 


> 


computer  file  structure  designed,  and  software  developed  to  allow 
data  from  digital  tapes  and  reports  to  be  loaded  into  the  base. 
This  effort  is  described  further  in  another  report^  and 
summarized  in  Section  2. 

Next,  data  availability  has  to  be  characterized  and  the 

observations  grouped  into  provinces  that  reduce  the  uncertainity 

of  the  physical  inputs.  Our  plan  is  to  group  the  data  a  priori 

2 

by  sediment  type  following  Hamilton's  classification  schemes. 
Since  we  concentrate  on  shallow  water  and  continental  margin 
areas,  however,  all  our  data  fall  into  one  physiographic  province 
under  Hamilton's  scheme  —  continental  terrace.  From  that  point 
the  only  classification  parameter  remaining  is  sediment  size 
class  (sand,  silt,  clay,  and  mixtures).  We  further  separate  our 
observations  based  on  an  analyses  of  variance  of  one  of  the 
physical  parameters  in  our  data  base  --  void  ratio  (an  expression 
of  porosity).  We  choose  void  ratio/porosity  because  other 
work9,10’16’17  shows  the  sensitivity  of  normal  incidence 
reflection  to  this  parameter.  After  provincing,  the  variability 
of  PHYSED  inputs  must  be  quantified  by  province.  This  step  is 
described  in  Section  2  and  in  reference  18. 

__ 

Molinelli,  E.J.,  D.A.  Waugh,  O.P.  Council,  and  B.A.  Brunson, 
1985,  "PHYPROSE  —  A  Data  Base  for  Physical  Properties  of 
Ocean  Sediments,"  PSI  TR-291304  for  NORDA  Code  113,  Planning 
Systems  Inc.,  McLean,  VA  22102. 

Molinelli,  E.J.,  B.A.  Brunson,  and  D.A.  Waugh,  1984,  "Acoustic 
Reflectivity  and  Mine  Burial  Properties  of  Sediments",  PSI  TR- 
301298  for  NCSC  Code  003,  Planning  Systems  Inc.,  McLean,  VA 
22101. 

1 7 

Brunson,  B.A.,  E.G.  McLeroy,  C.W.  Holland,  and  R.K.  Hagg,  1985, 
"Acoustic  Sea  Bottom  Classification:  A  Requirements 
Analysis",  PSI  TR-335313  for  NCSC  Code  401,  Planning  Systems 
Inc.,  McLean,  VA  22101. 

*  Council,  O.P. ,  R.K.  Hagg,  D.A.  Waugh,  and  E.J.  Molinelli,  1985, 
"Statistical  Analysis  of  the  PHYPROSE  Data  Base",  PSI  TR- 
291315  for  NORDA  Code  113,  Planning  Systems  Inc.,  McLean,  VA 
22101. 
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The  variability  of  the  physical  properties  within  a  province 
must  then  be  propagated  through  the  PHYSED  model  to  produce 
variability  in  the  geoacoustic  profiles  output.  We  choose  not  to 
attempt  to  propagate  errors  analytically  for  several  reasons. 
Generating  the  analytic  partial  derivatives  of  each  output  with 
respect  to  each  input  is  a  tedious  and  error  prone  operation, 
especially  because  of  the  use  of  complex  numbers;  the  resulting 
partials  will  usually  be  functions  of  the  other  variable  inputs 
and  sometimes  will  be  nonlinear  functions;  and,  there  is  an 
alternative  approach.  The  alternative  we  use  is  to  simulate  the 
variance  in  the  geoacoustic  profiles  by  randomly  sampling  from 
the  distribution  of  physical  properties  defined  for  a  province, 
performing  the  PHYSED  calculations  on  each  sample,  and 

accumulating  the  resulting  geoacoustic  profiles.  This  we  call 
the  Monte  Carlo  approach. 

For  the  Monte  Carlo  approach  several  physical  properties  are 
allowed  to  vary  independently  while  others  are  held  fixed.  We 
specifically  fixed  the  three  pore  fluid  properties  so  that  our 
variations  could  be  attributed  to  sediment  properties  and  not 
water  column  variations.  Also  fixed  were  the  structure  factor 
and  the  bulk  modulus  of  the  grains  for  which  no  data  were 

available  and  which  have  little  impact  on  compress ional  speeds. 

Void  ratio  (porosity),  grain  density,  and  mean  specific  surface 
(which,  with  porosity,  determines  permeability  and  the  pore  size 
parameter)  were  sampled  from  probability  distribution  functions 
generated  from  the  data  base.  Poisson's  ratio  (which  affects 

frame  bulk  modulus),  and  frame  log  decrements  (which  affect  frame 
losses)  were  varied  uniformly  over  a  conservative,  large  range 

O  1  Q 

based  on  values  published  by  Hamilton  ’  because  they  were  not 

19 

Hamilton,  E.L.,  1976,  "Attenuation  of  Shear  Waves  in  Marine 
Sediments",  J.  Acoust.  Soc .  Am.  60 ,  pp  334-338. 


available  in  the  data  base,  yet  had  shown  some  effect  on  BL  in 
previous  work.®  Frame  shear  modulus  was  derived  from  porosity 
using  the  "Stoll  stress"  procedure.1®  Where  approximations  had 
to  be  made  they  were  chosen  conservatively,  so  as  to  reduce  the 
variability  of  shallow  sediments  overall  or  to  increase  the 

variability  within  a  province.  In  this  way,  any  advantages 
discovered  for  the  provincing  will  be  robust,  i.e.  not  likely  to 
dissipate  under  different  approximations. 

A  further  approximation  is  the  use  of  "Stoll  stress" 

constants  appropriate  for  sands  for  the  case  of  shallow  sediments 
overall.  We  consider  this  simplification  preferable  to 

introducing  an  arbitrary  discontinuity  of  frame  modulus  to 

shallow  sediments  that  are  otherwise  continuous  in  void  ratio  and 
grain  size.  This  step  is  described  in  Section  3. 

The  PHYSED  model  must  be  run  at  single  frequencies.  Because 
the  frequencies  of  interest  are  those  associated  with  ASW 
detection,  we  select  5  frequencies  for  the  4  octaves  in  the 
range  100  to  1600  Hz. 

The  next  step  is  to  compute  bottom  loss  for  each  geoacoustic 
profile  resulting  from  the  PHYSED  Monte  Carlo  runs.  We  use  the 
numerical  model  REFLEC  to  perform  this  computation.  This  model 
has  been  developed  by  scientists  at  the  Naval  Ocean  Research  and 
Development  Activity  using  an  approach  described  by 
Brekhovskikh2®  and  has  been  used  in  studies  designed  to  test  the 
sensitivity  of  the  complex  reflection  coefficient  to  sediment 

O  I 

layering.  The  model  is  capable  of  using  the  compressional  wave 

20 

Brekhovskikh,  L.M. ,  1960,  Waves  in  Layered  Media,  Academic 
Press,  New  York. 

^Gilbert,  K.E.,  1980,  "Reflection  of  Sound  from  a  Randomly 
Layered  Ocean  Bottom",  J.  Acoust.  Soc .  Am..  68 ,  pp  1454-1458. 
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properties  provided  by  the  PHYSED  model  and  generating  a  multi¬ 
layered  approximation  to  the  vertical  profiles.  The  number  and 
thickness  of  the  layers  are  selected  automatically  by  a 
preprocessor  which  takes  into  account  the  sound  speed  gradient 
and  the  frequency  at  which  the  reflection  coefficient  calculation 
is  to  be  performed.  The  layers  are  not  allowed  to  exceed  0.1 
wavelengths  in  thickness,  thus  ensuring  a  smooth  approximation  of 
the  sound-speed  profile.  The  output  of  REFLEC  consists  of  either 
the  complex  reflection  coefficient  or  a  bottom  loss.  These  are 
available  at  any  grazing  angle  specified  for  any  frequencies  of 
interest.  The  model  can  produce  estimates  for  any  desired 
frequency  bandwidth  (e.g.,  one-third  octave)  by  averaging  the 
results  for  discrete  frequencies  within  the  desired  band; 
however,  we  do  not  employ  the  bandwidth  averaging  option.  REFLEC 
does  not  account  for  energy  loss  due  to  the  generation  of  shear 
waves,  hence  REFLEC  is  a  "fluid"  sediment  model  not  a  "solid" 
sediment  model.  We  find  REFLEC  satisfactory  for  our  purposes  as 
discussed  in  Section  4.  For  the  problem  of  long-range  detection, 
the  BL  at  the  low  grazing  angles  are  of  primary  interest.  We 
decide  that  a  grazing  angle  of  5°  is  representative  of  BL  that 
impacts  detection.  The  generation  of  BL  distributions  is 
described  in  Section  4. 

Finally,  we  assess  our  results  in  terms  described  in  the 
objective  --  mean  and  variance  of  BL  in  an  individual  province 
and  for  shallow  sediments  overall,  difference  between  an 
individual  province  and  all  shallow  sediments,  accuracy  of 
predictions,  and  availability  of  information  to  assign  a 
geographic  position  to  a  province.  This  step  requires  careful 
statistical  procedures  and  the  analyses  of  variance  as  described 
in  Section  5.0. 


There  are  a  few  assumptions  built  into  our  approach.  We 
cannot  know  exactly  the  distribution  of  a  property  in  a  province; 


we  can  only  estimate  the  distribution  based  on  samples  from  that 
province.  We  therefore  use  statistical  methods  to  help  us 
distinguish  those  conclusions  which  can  be  drawn  with  a  high 
probability  of  being  correct.  Stated  another  way,  our  values  for 
means,  variances,  differences,  probabilities,  and  percentages  are 
estimates  in  which  we  have  some  statistical  confidence;  but  these 
values  can  not  be  considered  exact. 

We  also  make  the  assumption  that  the  marine  sediment 
physical  data  we  could  assemble  into  a  data  base  with  one  man- 
year  of  effort  represents  the  sediment  physical  data  available  at 
large.  We  have  made  every  effort  to  obtain  a  cross  section  of 
physical  observations  on  the  continental  terrace  (continental 
shelf,  continental  slope,  and  continental  rise)  so  we  further 
assume  the  data  we  did  assemble  is  representative  of  continental 
terrace  conditions.  Though  the  conclusions  we  reach  can  not  be 
applied  to  continental  margins  everywhere,  we  believe  that  to 
justify  the  large  effort  required  to  do  that  complete  job,  it  is 
first  necessary  that  the  data  set  we  have  here  assembled  show 
some  exploitable  trends. 


In  our  calculation  of  the  distribution  of  properties  from 
the  data  base  we  assume  measurements  from  different  depths  in  the 
same  core  are  independent  observations.  Each  depth  is 
independently  assigned  to  a  province.  This  can  skew  our 
distributions  because  deep  cores  come  from  soft  sediments  and 
thus  more  weight  is  assigned  per  soft  core  than  to  a  hard 
sediment  core.  When  we  generate  a  geoacoustic  profile  during  a 
Monte  Carlo  run,  the  only  depth  dependence  is  that  associated 
with  the  Stoll  stress  procedure;  porosity  and  the  other 
properties  identified  for  the  run  remain  constant  with  depth.  We 
thus  are  modeling  a  vertically  homogeneous  sediment  that  is  one 
member  of  the  class  of  sediments  belonging  to  a  province. 


2.0  THE  PHYSICAL  DATA 


2.1  The  PHYPBOSE  Data  Base 

As  described  in  a  separate  report, all  the  physical  sedi¬ 
ment  data  used  in  this  study  were  obtained  from  the  National  Geo¬ 
physical  Data  Center  (NGDC)  in  Boulder,  Colorado.  NGDC  supports 
a  computerized  summary  of  their  digital  and  technical  report 
holdings  which  was  searched  for  us.  We  eliminated  from  our 
search  those  holdings  that  had  only  descriptive  accounts  of 
sediment  color  or  biological  stratification.  The  physical 
properties  upon  which  we  could  key  our  searches  were  sediment 
"texture"  data,  "engineering"  data  and  "acoustic"  data.  By  far 
the  most  extensive  holdings  (thousands  of  entries)  were  found  for 
texture  data.  The  descriptions  of  the  texture  holdings  were 
studied  to  identify  those  which  included  "raw"  texture  data  and 
simultaneously  had  acoustic  or  engineering  data.  We  noted 
whether  measurements  were  in  high-priority,  shallow-water,  areas. 
When  we  had  thus  narrowed  our  choices,  we  ordered  data  to  cover 
as  many  areas  as  possible,  favoring  digital  data  sets  with  large 
numbers  of  stations  over  similarly  located  data  published  in 
reports.  Digital  data  were  ordered  as  magnetic  computer  tape; 
non-digital  data  were  ordered  as  microfiche. 

Using  the  criteria  set  forth  above,  twenty-seven  (27)  data 
reports  were  ordered  from  NGDC  comprising  95  sheets  of  micro¬ 
fiche.  The  area  represented  by  each  data  report,  its  NGDC  iden¬ 
tification  number,  the  number  of  microfiche  sheets,  remarks  on 
its  contents  and  disposition  with  respect  to  the  PHYPROSE  data 
base,  are  summarized  in  Table  2-1. 

The  data  entry  process  consisted  of  reading  the  microfiche 
and  tabulating  PHYPROSE  parameters  onto  work  sheets.  Then  a 
BASIC  language  program  on  a  TEKTRONIX  4052  desktop  computer  was 
run  to  accept  information  in  many  different  units  and  combina¬ 
tions  and  to  output  formatted  records  using  PHYPROSE  conventions. 
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NGDC 

NO.  Of 

Area 

£0  No. 

Sheets 

0.S.  Meat  Coaat 

09005007 

6 

V, 

09005058 

5 

•- 

09595004 

2 

£ 

09825001 

2 

20995001 

3 

*  " 

rv. 

Mediterranean  Sea 

09005022 

1 

09055005 

6 

09535002 

2 

i 

03045003 

3 

Arctic 

09245001 

3 

09325001 

1 

K 

09375002 

6 

c-: 

09995032 

4 

09995004 

2 

b 

Persian  Sea 

09315004 

8 

Car lbbean 

03035001 

4 

09005012 

3 

09065007 

3 

09265005 

2 

s' 

09295001 

2 

Gulf  of  Alaska 

06255002 

2 

*r 

Indonesia 

09255001 

1 

Indian  Ocean 

09005024 

1 

% 

09645001 

5 

09785006 

5 

,s 

s 

09785007 

7 

So.  Hemisphere 

09995010 

5 

a 


NO.  of 
Stations 

Entered 

Ra narks 

39 

Texture,  Mineralogy, 

6  Engineering 

102 

X 

NCEL,  Texture, 

6  Engineering 

10 

X 

NOO,  Texture,  Engin¬ 
eering,  C  Acoustics 

10 

NOO,  Texture,  Engin¬ 
eering,  6  Acoustics 

1.103 

Texture 

3 

NOO,  Texture,  6 
Engineering 

21 

X 

NOO,  Texture,  Engin¬ 
eering,  6  Acoustics 

59 

NOO,  Texture 

7 

X 

NOAA,  Texture,  6 
Engineering 

64 

NOO,  Texture 

38 

USSR 

54 

NOO,  Texture  6 
Engineering 

20 

X 

NOO,  Texture  t 
Engineering 

154 

X 

NOO,  Texture 

108 

NOO,  Texture 

31 

X 

Engineering 

65 

X 

Texture,  Engineering 

26 

NOO,  Texture  6 
Engineering 

12 

X 

NOO,  Texture,  Engin¬ 
eering,  6  Acoustics 

40 

NOO,  Texture 

16 

X 

Texture,  Engineering 

27 

NOO,  Texture 

8 

NOO,  Texture 

57 

X 

NOO,  Texture 

4 

X 

NOO,  Texture,  Engin¬ 
eering,  6  Acoustics 

9 

Too  dark  to  read 

33 

No  data 
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These  records  were  transferred  to  a  VAX  11/780  run  by  the  Naval 
Ocean  Research  and  Development  Activity  (NORDA)  at  NSTL  Station, 
Mississippi  through  a  a  communications  bus  and  telephone  connec¬ 
tion.  Listings  of  the  resultant  VAX  files  were  made  and  compared 
to  the  original  microfiche  for  quality  control.  Finally,  VAX 
FORTRAN  programs  were  run  to  load  the  records  into  PHYPROSE. 
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In  Table  2-1,  no  mark  in  the  "Entered"  column  indicates  that 
data  from  that  report  have  not  been  entered  into  PHYRPOSE.  Only 
eight  of  twenty-seven  reports  have  been  incorporated  because  the 
data  must  be  entered  manually  and  the  process  requires  extensive 
uses  of  manpower.  The  reports  were  added  to  the  data  base  in  an 
order  that  reflected  the  importance  of  the  area,  the  widest  geo¬ 
graphic  coverage  possible,  and  the  expectation  of  additions  from 
digital  sources.  Consequently  no  U.S.  east  coast  data  were 
entered  by  hand  (or  even  ordered)  because  of  the  vast  number  of 
stations  available  from  the  USGS  tape  discussed  below.  Substan¬ 
tial  numbers  of  stations  were  expected  from  digital  tapes  for  the 
Arctic,  the  Gulf  of  Mexico,  the  Pacific  margin,  and  southeast 
Asia.  Also  anticipated  was  a  digital  tape  of  unclassified 
NAV0CEAN0  data;  so,  microfiche  of  NAVOCEANO  observations  were  not 
used  until  the  contents  of  that  tape  could  be  viewed  to  avoid 
duplication  of  effort.  When  the  tape  proved  extremely  limited  in 
coverage  (due  to  classification)  late  in  the  effort  only  a  few 
NAVOCEANO  stations  could  be  entered. 


4  , 


•  .  "  v 

•  • w  * 


During  the  search  of  NGDC  holdings,  several  digital  computer 
tapes  were  identified  as  having  useful  data  in  priority  areas. 
Five  of  these  tapes  were  ordered  to  provide  as  wide  a  global 
coverage  as  possible.  The  area  represented  by  each  tape,  its 
NGDC  identification  number,  the  number  of  stations  it  contains, 
its  disposition  with  respect  to  PHYPROSE,  and  remarks  on  its  con¬ 
tent  are  presented  in  Table  2-2.  The  NAVOCEANO  tape  was  not 
ordered  from  NGDC  but  received  directly,  therefore  it  is  not 
included  in  this  table.  Because  it  did  not  prove  useful,  it  will 
not  be  discussed  further. 
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Table  2-2.  Computer  Tapes  Ordered  from  NGDC 


AREA 

NGDC  ID  # 

STATIONS 

ENTERED 

REMARKS 

U.S.  East  Coast 

06995002 

3715 

yes 

USGS,  texture 

Arctic 

20995002 

688 

yes 

U.  Washington, 
texture 

Gulf  of  Mexico 

23055001 

561 

yes 

Texas  A&M, 
texture 

North  Pacific 

15995012 

688 

no 

N00( Scripps ) -- 
no  useful  data 

Southeast  Asia 

15995014 

8168 

no 

N00( Scripps ) -- 
no  useful  data 

The  process  of  retrieving  data  from  the  computer  tapes  con¬ 
sisted  of  writing  and  running  special  purpose  VAX  FORTRAN  soft¬ 
ware  to  read  data  from  the  tapes  and  load  the  information  into 
PHYPROSE-like  files.  These  were  screened  for  errors,  edited  and 
appended  to  the  end  of  the  actual  PHYPROSE  files.  Not  all  sta¬ 
tions  on  the  tape  were  appended  to  PHYPROSE  because  not  all  sta¬ 
tions  had  values  for  the  select  PHYPROSE  parameters.  Especially 
disappointing  was  the  fact  that  the  NAVOCEANO  (N00)  data  compiled 
onto  tape  by  the  Scripps  Institution  of  Oceanography  (the  last 
two  tapes  in  Table  5-2)  contained  no  information  other  than  qual¬ 
itative  grain  size. 

2.2  PHYPROSE  Contents 

The  contents  of  the  PHYPROSE  data  base  as  of  1  April  1985 
are  described  in  terms  of  the  total  number  of  stations,  the 
number  of  stations  with  particular  types  of  data  (common  proper¬ 
ties,  engineering  properties,  etc.),  the  total  number  of  records 
of  each  data  type  and  maps  of  station  locations.1'’  In  summary, 
PHYPROSE  contains: 
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*4,355  stations 

*  233  stations  reporting  common  properties  (including  void 
ratio  &  grain  density) 

2,277  records 

*  106  stations  reporting  engineering  properties  (shear 
strength,  etc.) 

933  records 

*1,950  stations  reporting  mineral  properties  (percent 
Calcium  Carbonate,  etc.) 

2,935  records 

*2,109  stations  reporting  qualitative  grain  sizes  (sand, 
silty  sand,  etc.) 

2,898  records 

*  180  stations  reporting  size  class  information  (percent 
sand,  percent  silt,  etc.) 

1,353  records 

*3,900  stations  reporting  full  grain  size  distributions 
(percent  by  weight  in  one  phi  size  bins) 

9,214  records 

*  40  stations  reporting  acoustic  properties  (including 
compressional  speed) 

921  records 

*  40  stations  reporting  fluid  properties  (including 
temperature)  168  records 

These  stations  are  distributed  over  17  ocean  areas  as 
follows  (ocean  area  code  in  bold  numerals): 

*U.S.  East  Coast  11  (Nova  Scotia  to  Key  West) 

2,109  stations 

♦Caribbean  Sea  12  (Lesser  Antilles) 

48  stations 

♦Gulf  of  Mexico  13  (Mississippi  delta) 

206  stations 
♦Beaufort  Sea  24 


13  stations 


♦Chukchi  Sea  25 

700  stations 
♦East  Siberian  Sea  26 
154  stations 
♦Laptev  Sea  27 

114  stations 
♦Barents  Sea  29 

2  stations 

♦Mediterranean  Sea  30,  31 
104  stations 
♦Black  Sea  34 

2  stations 

♦Indian  Ocean  margin  50,  51 
13  stations 

♦U.S.  West  Coast  60  (Los  Angeles  &  Gulf  of  Alaska) 

101  stations 
♦South  China  Sea  62 
20  stations 
♦Sea  of  Okhotsk  66 
20  stations 
♦Bering  Sea  shelf  67 
739  stations 

One  of  the  most  dramatic  features  of  the  PHYPR0SE  holdings 
is  the  vast  quantities  of  textural  (i.e.,  grain  size)  information 
and  the  lack  of  geoacoustic  data.  It  is  true  that  a  bias  against 
acoustic  data  resulted  from  our  decision  to  pass  over  NAV0CEAN0 
reports  ("lab  items")  in  anticipation  of  receiving  a  digital 
tape.  However,  had  we  biased  our  data  base  in  the  opposite 
direction  by  incorporating  all  stations  with  acoustic  data  from 
our  twenty-seven  reports,  we  would  have  had  only  57  stations  with 
geoacoustic  measurements  --  only  1%  of  the  stations  that  are  in 
the  data  base  already,  and  less  than  1%  of  the  stations  had  we 
entered  all  the  data  listed  in  Table  2-1.  It  is  fair  to  state 
that  there  is  a  vast  amount  of  textural  data  and  a  lack  of 
geoacoustic  data.  However,  for  PHYPR0SE  to  enable  testing  of 
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Biot/Stoll  outputs  some  geoacoustic  data  had  to  be  included  in 
the  data  base. 

The  ratio  of  the  number  of  texture  observations  to  the 
number  of  common  properties  observations  (including  void  ratio 
and  grain  density)  also  seems  well  represented  by  PHYPROSE.  In 
the  listing  of  high  quality  microfiche  in  Table  2-1,  there  are 
about  2040  stations  described;  2009  of  them  report  texture  data. 
In  the  same  set,  about  420  stations  report  engineering  properties 
(including  common  properties).  Therefore,  in  the  specialized 
subset  described  in  Table  2-1,  common  properties  are  available 
about  20%  of  the  time  and  texture  properties  are  available  98%  of 
the  time.  Given  that  Table  2-1  selectively  ignored  purely  tex¬ 
tural  data  sets,  one  can  expect  common  properties  to  be  available 
much  less  than  20%  of  the  time.  In  PHYPROSE,  common  properties 
are  present  5%  of  the  time.  This  is  not  unreasonable. 
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We  therefore  conclude  that  PHYPROSE  adequately  represents 
the  relative  abundance  of  geoacoustic  data,  void  ratio  and  grain 
density  data,  and  grain  size  distribution  data  in  shallow  waters. 
The  stated  abundances  indicate  that  a  great  advantage  in  predict¬ 
ing  BL  is  achieved  if  a  geographic  position  in  shallow  water  can 
be  allocated  to  a  BL  province  based  on  grain  size  data  rather 
than  on  acoustic  data,  geoacoustic  data,  or  engineering  data. 

2.3  The  Definition  of  Provinces 

It  is  beyond  the  scope  of  this  effort  to  discover  the 
optimum  provincing  of  ocean  sediments  based  on  BL.  We  wish  first 
to  demonstrate  that  there  is  some  advantage  to  provincing;  speci¬ 
fically  we  wish  to  quantify  the  advantage  in  provincing  based  on 

O 

the  Hamilton  classification  scheme.  For  shallow  water  sediments 
on  the  continental  shelf,  continental  slope,  and  continental 
rise,  Hamilton  defines  a  single  physiographic  province  called  the 
continental  terrace.  Within  the  terrace,  he  separates  sediments 
into  qualitative  size  classes  --  coarse  sand,  fine  sand,  very 
fine  sand,  silty  sand,  sandy  silt,  silt,  sand-silt-clay,  clayey 
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silt,  and  silty  clay.  Based  on  212  samples  assigned  to  these 
size  classes  he  describes  ranges  for  physical  and  geoacoustic 
parameters.  Hamilton's  scheme  is  shown  in  Figure  2-1. 

For  the  large  amount  of  data  in  PHYPROSE  we  need  to  define  a 
qualitative  sediment  size  class  based  on  grain  size  data  unambig- 
iously.  By  convention,  the  gravel  size  class  includes  grains 
with  diameters  greater  than  2  millimeters.  A  convenient  scale 
for  grain  diameters  is  the  phi-scale.  By  definition: 

<J>  =  -logg  (grain  size  in  millimeters) 

On  the  phi  scale,  gravel  corresponds  to  values  less  than  -1.  The 
sand  size  class  includes  grains  with  diameters  from  2  down  to 
.0625  millimeters(-l<  <j>  <4);  and,  silt  grains  are  in  the  range 
from  .0625  to  .0039  millimeters  (4<<t><8).  The  clay  size  class 
includes  all  grains  smaller  than  .0039  millimeters  (<t>>8).  The 
relative  percentages  of  the  sediment  in  these  size  classes  can  be 
used  to  define  an  unambiguous,  qualitative  grain  size  class  as 
illustrated  in  Figure  2-2.  We  adopt  a  numbering  system  for  the 
qualitative  size  classes  as  shown  in  that  figure;  sand  is  sedi¬ 
ment  size  class  1,  silty  sand  is  2,  etc.  Not  shown  is  an 
eleventh  class,  gravel,  which  is  presumed  to  describe  sediments 
with  more  than  50%  gravel  fraction.  No  sediments  in  the  PHYPROSE 
data  base  were  so  labeled.  The  percentages  that  bound  these 
classes  are  consistent  with  the  data  displayed  by  Hamilton. ^ 
Silty-clay  mixtures  are  termed  muds  in  our  convention;  and,  no 
distinction  is  made  among  the  sands  (coarse,  fine,  very  fine)  in 
our  convention. 

The  separation  of  continental  terrace  sediments  by  qualita¬ 
tive  size  class  is  an  a  priori  provincing  likely  to  affect  BL 
because  of  the  differing  compress ional  speeds  claimed  for  each 
type  in  Figure  2-1.  Our  experience  with  PHYSED  modeling  has 
shown  the  importance  also  of  poros ity/void  ratio  to  bottom 
reflectivity . ^  We  therefore  tested  whether  continental 
terrace  sediments  could  be  separated  based  on  their  void  ratio 


TABLE  IB.  Continental  terrace  (shelf  and  slope)  environment;  sodlmoat  densities,  porosities, 
sound  velocities,  sod  velocity  ratios. 


Density* 

Porosity  * 

Velocity* 

fr/enr) 

l*> 

Un/s) 

V«  toe  tty  ruto* 

Sod  lm  eat  type 

Av  SE 

At  SE 

At  SE 

At  SE 

Seed 


Coars# 

2.034 

38.8 

... 

1838 

... 

1.201 

Floe 

1.941 

0.023 

45.6 

1.02 

1748 

11 

1.145 

0.006 

Very  fine 

IJS6 

0.022 

$0.0 

0J7 

1702 

18 

1J1S 

0.012 

Silty  sand 

1.722 

0.020 

$53 

0.72 

1648 

10 

1.078 

0.006 

Sandy  silt 

1.771 

0.033 

$4.1 

1.49 

1652 

12 

1.080 

0.007 

SiU 

1.740 

0.047 

S6.3 

1  JO 

1615 

a 

1.057 

0.005 

San  d-e  tit-clay 

IJ96 

0.022 

8«  J 

1.53 

1579 

8 

1.033 

0.005 

Clayey  silt 

1.488 

0.016 

71.6 

9J6 

1S49 

4 

1.014 

0.303 

Silty  clay 

1*421 

0.01S 

7S.9 

0.82 

1520 

3 

0.994 

0.002 

‘Laboratory  values:  23  U,  1  atm;  density:  Saturated  bulk  density;  porosity:  Salt-free;  velocity 
ratio:  Velocity  in  sediment/ velocity  uj  sea  water  at  23  *C,  l  atm,  and  salinity  of  sediment  pore 
water.  SE:  Standard  error  of  the  eieae. 


Figure  2-1.  From  Hamilton  ,  page  1315 


V 


Figure  2-2.  The  classification  of  sediment  by  relative 

\.J  proportions  of  grain  sizes.  A  sediment  that  is  more  than  75% 

sand  is  classed  as  sand;  a  sediment  that  is  less  than  25"  sand  is 
.v  classed  as  silt,  mud  or  clay  depending  upon  the  non-sand  portion. 

^  If  the  non-sand  portion  is  more  than  66%  silt,  the  sediment  is 

termed  a  silt;  between  66%  and  33%  silt,  the  sediment  is  a  mud; 
,  and,  less  than  33%  silt,  the  sediment  is  a  clay.  Mixtures  are 

classified  according  to  the  diagram. 


values,  and  found  that  both  depth  of  the  seafloor  and  ocean  area 
were  additional  criteria  more  important  than  sediment  size  type 
in  explaining  void  ratio  variance. 1®  We  summarize  the  results  of 
that  study  in  the  next  paragraphs. 

Figure  2-3  shows  a  histogram  of  all  void  ratio  values  in  the 
PHYPROSE  data  base  for  which  there  was  also  grain  size  infor¬ 
mation.  The  distribution  shows  two  major  peaks  and  is 
asymmetric,  strongly  suggesting  that  there  are  several  factors 
influencing  void  ratio.  We  use  one-way  classification  analysis 
of  variance  (ANOVA)  to  demonstrate  the  relative  significance  of 
four  "treatments",  i.e.  sediment  grain  size  type,  seafloor  depth, 
ocean  area,  and  depth  in  sediment  core.  For  the  sediment  type 
"treatment"  the  ten  classes  already  described  and  depicted  in 
Figure  2-2  are  used.  For  the  seafloor  depth  treatment  we  define 
four  classes:  shelf  (0  to  200  m),  upper  slope  (200  to  1000  m), 
lower  slope  (1000  to  2000  m) ,  and  continental  rise  (2000  to  4000 
m).  For  ocean  area  we  use  the  regions  defined  for  the  PHYPROSE 
data  base;  the  ocean  areas  with  data  were  listed  in  Section  2.1. 
Finally,  we  define  depth  intervals  within  the  core  as  follows:  0 
to  10  cm,  10  to  60  cm,  60  to  80  cm,  80  to  100  cm,  100  to  120  cm, 
120  to  140  cm,  and  deeper  than  140  cm. 

One-way  ANOVA  calculations  for  each  treatment  are  presented 
in  Tables  2-3  through  2-6.  Here,  the  influence  of  water  depth, 
sediment  type,  ocean  area  and  depth  in  sediment  core  on  void 
ratio  are  separately  tested.  The  table  shows  the  F-statistic  for 
each  factor.  These  values  are  compared  to  the  distribution  of 
the  F-statistic  for  the  appropriate  degrees  of  freedom  and  all 
four  values  are  found  to  be  significant.  This  means  that  there 
is  a  significant  difference  in  void  ratio  between  at  least  one 
pair  of  treatment  classes  in  each  of  the  four  treatments;  it  does 
not  guarantee  that  all  treatment  classes  are  significantly  dif¬ 
ferent  from  all  others. 
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Figure  2-3.  Histogram  of  all  void  ratio  values 


Table  2-3.  One-way  ANOVA  on  Seafloor  Depth 


Depth 

Number  of 
Observations 
(n) 

Sum 

Sum  of 
Squares 

Contribution 
to  Square  of 
Sum  per 
Observation 

Mean 

0-200 

237 

324.021 

480.1941 

442.994 

1.367 

200-1000 

227 

556.130 

1571.374 

1362.469 

2.449 

1000-2000 

136 

429.288 

1634.992 

1325.060 

3.156 

2000 

492 

1870.780 

7879.046 

7113.451 

3.802 

=  337.419 


=  1.187 


F  =  284.224 


Table  2-4.  One-way  ANOVA  on  Ocean  Area 


103.479 


1 . 36167 


75.99 
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Table  2-5.  One-way  ANOVA  on  Depth  in  Sediment  Core 


Sediment 

Core 

Depth  (m) 

Number  of 
Observations 
(n) 

Sum 

Sum  of 
Squares 

Contribution 
to  Square  of 
Sum  per 
Observation 

Mean 

0-.10 

59 

115.515 

296.861 

226.165 

1.957 

0.10-0.60 

215 

455.493 

1199.897 

964.995 

2.118 

0.60-0.80 

59 

123.004 

331.719 

256.440 

2.085 

0.80-1.00 

27 

57.399 

154.258 

122.024 

2.126 

1.00-1.200 

6 

14.700 

45.442 

36.015 

2.450 

1.200-1.400 

9 

23.566 

80.403 

61.706 

2.618 

1.400 

17 

55.832 

192.016 

183.365 

3.284 

Sc 

2  =  4.5038 

S2  = 

1.1685 

F  =  2.85 

Table  2-6. 

One-way  ANOVA  on  Sediment  Type 

Sediment 

Number  of 
Observations 

Sum  of 

Contribution 
to  Square  of 
Sum  per 

Type 

(n) 

Sum 

Squares 

Observation 

Mean 

Sand 

12 

11.66 

11.605 

11.337 

0.972 

Si  lty 

Sand 

80 

223.783 

855.108 

625.734 

2.797 

Muddy 

Sand 

51 

118.007 

350.074 

273.052 

2.3139 

Sandy 

Silt 

57 

81.225 

135.150 

115.746 

1.425 

Sandy 

Mud 

48 

78.0110 

140.938 

126.786 

1.625 

Sandy 

Clay 

4 

7.11 

13.153 

12.638 

1.777 

Silt 

87 

172.372 

396.047 

341.518 

1.981 

Mud 

314 

906.089 

3220.703 

2614.64 

2.886 

Clay 

439 

1581.823 

6442.821 

5699.69 

3.603 

S  2  =  62.7677  S2  = 
c 


F  =  38.96 
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Inspection  of  Tables  2-3  through  2-6  indicates  that  the  mean 
square  between  classes,  Sc^>  is  largest  for  variation  due  to 
water  depth.  This  suggests  that  water  column  depth  has  a  greater 
influence  on  void  ratio  than  the  other  factors.  The  influence  of 
water  depth  on  void  ratio  was  further  examined  by  plotting  histo¬ 
grams  of  void  ratio  for  4  different  water  depth  intervals.  These 
plots  are  presented  in  Figure  2-4.  The  distributions  at  1  to  200 
and  greater  than  2000  m  seem  reasonably  homogeneous  (although 
this  was  not  explicitly  tested).  The  distributions  for  the  con¬ 
tinental  slope  (200-1000  m  and  1000-2000  m)  in  contrast  show  two 
separate  peaks,  i.e.  are  bimodal,  suggesting  that  at  these  depths 
some  other  factor  influences  void  ratio.  These  two  distributions 
are  therefore  divided  up  according  to  ocean  area  code  --  that  is, 
according  to  the  area  in  the  ocean  where  the  measurements  were 
taken.  The  void  ratios  in  the  200-1000  m  interval  come  from  the 
Mediterranean  Sea  and  the  Gulf  of  Alaska  as  shown  in  Figure  2-5a 
and  2-5b.  The  void  ratios  in  the  1000  to  2000  m  interval  come 
from  the  Caribbean,  the  Mediterranean  and  the  South  China  Seas  as 
shown  in  Figures  5c,  d  and  e.  The  high  void  ratios  in  the  deeper 
interval  are  associated  with  the  void  ratios  for  the  South  China 
Sea  shown  in  Figure  2-5e.  Also,  the  high  void  ratio  data  in  Fig¬ 
ure  2-4b  for  the  200-1000  m  interval  seems  to  come  primarily  from 
void  ratio  taken  in  the  Gulf  of  Alaska,  shown  in  Figure  2-5b. 

These  results,  however,  still  support  the  general  interpre¬ 
tation  that  depth  has  a  greater  influence  on  void  ratio  than 
ocean  area.  The  separation  in  the  void  ratio  data  created  by 
grouping  according  to  seafloor  depth  is  greater  than  that  created 
by  grouping  according  to  ocean  area  in  Figure  2-5.  More  gener¬ 
ally,  the  histograms  in  Figure  2-4  indicate  that  no  other  factor 
is  likely  to  have  as  great  an  influence  on  the  void  ratio  data  as 
water  depth.  A  two-way  classification  ANOVA  calculated  for  the 
effects  of  sediment  type  and  seafloor  depth  simultaneously  gives 
a  similar  result.  Seafloor  depth  has  a  greater  influence  on  void 
ratio  than  sediment  grain  size  class. 


water  column  depth  l-200m 


b  -  water  column  depth  200-1000m 


ora  in  in  X7»  in  in  m  47a 


an  iTa  in  17a  *n 


water  column  depth  1000-2000m 


d  -  water  column  depth  2000m 


Figure  2-4.  Histogram  of  void  ratio  values  separated  by  water  column  depth 
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-  water  depth  200-1000m,  ocean  area  b  -  water  depth  200-1000m,  ocean  area 
=  western  Mediterranean  Sea  =  Gulf  of  Alaska 


c  -  water  depth  1000-2000m  ocean  area  d  -  water  depth  1000-2000m,  ocean  area 
=  Caribbean  Sea  =  western  Mediterranean  Sea 
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e  -  water  depth  1 000-2000m,  ocean  area  =  South  China  Sea 

Figure  2-5.  Histograms  of  void  ratios  for  the  two  depth  intervals  200-1000m 
and  1000-2000m  separated  according  to  ocean  area 


Because  of  the  importance  of  seafloor  depth  for  void  ratio, 
and  the  importance  of  void  ratio/porosity  for  seafloor  reflectiv¬ 
ity,  we  have  provinced  continental  terrace  sediments  by  seafloor 
depth  in  addition  to  sediment  grain  size  type.  To  keep  the 
number  of  provinces  manageably  small,  and  the  number  of  samples 
per  province  meaningfully  large,  we  have  combined  some  classes 
using  for  guidance  the  mean  values  listed  in  the  one-way  ANOVA 
tables  (2-3  and  2-6)  and  the  histogram.  Thus  we  combine  the 
upper  and  lower  slope  depths  into  one  class  (slope  --  200  to 
2000  m)  ;  and  we  group  silty  sand  with  muddy  sand  (and  call  the 
group  dirty  sand),  we  group  sandy  silt  with  sandy  mud  and  with 
sandy  clay  (sandy  muds),  and  we  group  mud  and  clay  (muds). 

We  do  not  separate  continental  terrace  sediments  by  ocean 
area  at  first  in  order  to  keep  sample  sizes  large;  but,  we 
reserve  these  ocean  area  criterion  for  use  during  refinemen 
stages  if  we  do  not  obtain  sufficient  separation  of  BL  with  the 
original  provincing  scheme.  As  will  be  shown  in  Section  5,  there 
is  adequate  separation  of  BL  for  the  purpose  of  this  report  using 
the  original  provincing  scheme. 

We  adopt  a  two-letter  naming  convention  to  facilitate 
references  to  the  provinces.  The  first  letter  refers  to  the 
seafloor  depth  and  the  second  letter  refers  to  the  grain  size 
type  using  the  codes  in  Table  2-7. 

2.4  Province  Physical  Properties 

There  were  sufficient  data  in  the  PHYPROSE  data  base  to 
describe,  for  each  shelf  and  slope  province,  the  distribution  of 
three  of  the  PHYSED  input  parameters  --  void  ratio/porosity, 
grain  density,  and  mean  specific  surface  (which  drives  permeabil¬ 
ity  and  pore  size  paramter).  Table  2-8  lists  the  number  of 
observations  available  from  PHYPROSE  for  each  of  the  three  param¬ 
eters  in  each  of  eleven  provinces.  Province  ZZ  is  taken  to 
represent  continental  terrace  sediments  as  a  whole.  Measurements 
at  different  depths  within  the  same  core  were  included  as 
separate  observations. 
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Table  2-7.  Province  Name  Code 


CT  V  V 

>v 

£ 


'  V  V  S 


* 


r.y 


First  Letter 

Second  Letter 

m 

(seafloor 

depth) 

(grain  size  type) 

Letter 

Depth 

Letter 

Name(s) 

Class 

(s) 

Interval  (m) 

Number 

*  *• 

A 

0-200 

A 

Sand 

1 

B 

Silty  sand 

2 

B 

201-2000 

Muddy  sand 

3 

Clayey  sand 

4 

> 

C 

2001-4000 

C 

Sandy  silt 

5 

Sandy  mud 

6 

V 

Sandy  clay 

7 

Z 

All  depths 

D 

Silt 

8 

3 

E 

Mud 

9 

Clay 

10 

V 

Z 

All  size 

1 

thru 

10 

types 

4 


Table  2-8.  Number  of  Observations  that  Define 
Distributions  for  the  Provinces 


Void  Ratio: 

Grain  Density: 

Mean  Specific 
Surface : 

Province 

Number  of 

Number  of 

Number  of 

Code 

Observations 

Observations 

Observations 

AA 

8 

16 

1,870 

y 

AB 

24 

48 

822 

AC 

36 

37 

1,061 

A 

>. 

AD 

22 

22 

2,359 

AE 

41 

43 

1,140 

BA 

2 

18 

169 

• . 

BB 

10 

17 

111 

BC 

44 

49 

164 

.V 

BD 

4 

5 

42 

BE 

294 

310 

451 

ZZ 

1,457 

1,542 

9,214 

*4 

.-'V 

■V 

IV.; 


v  v 


Tv;?: 

*; 

v- 

ft: 


The  distributions  are  represented  as  cumulative  distribution 
functions  which  are  more  easily  computed  from  real  data  than  his¬ 
tograms  because  no  bins  have  to  be  defined,  and  which  are 
directly  used  by  the  Monte  Carlo  sampling  routines.  The  cumula¬ 
tive  distribution  functions  are  plotted  for  void  ratio  in  Fig¬ 
ure  2-6,  for  grain  density  in  Figure  2-7,  and  for  mean  specific 
surface  in  Figure  2-8. 

These  distributions  were  computed  by  sorting  according  to 
value  the  observations  of  each  parameter  in  each  province, 
assigning  each  sorted  observation  a  sequence  number  (from  1  to 
the  total  number  of  observations),  defining  the  "probability"  by 
dividing  each  sequence  number  by  the  total  number  of 
observations,  dropping  all  duplicate  points  but  the  one  with  the 
highest  sequence  number,  and  defining  a  zero  probability  point  by 
extrapolating  the  mean  slope  back  to  the  parameter  axis 
(providing  it  did  not  intersect  below  the  smallest  possible  value 
of  the  parameter  —  zero  for  void  ratio  and  mean 
specific  surface,  one  for  grain  density).  The  resulting 
probability  associated  with  each  value  on  the  parameter  axis 
represents  an  estimate  of  the  fractional  chance  that  the  given 
parameter  value  or  a  lesser  one  will  occur  in  that  province. 

The  void  ratios  in  Figure  2-6  have  about  a  value  of  1.0  in 
the  coarse  sediments  and  vary,  predominately  between  1.0  and  2.0 
in  the  finer  sediments,  except  for  muds  (Provinces  AE  and  BE) 
which  usually  have  values  above  2.0.  Terrace  sediments  overall 
(Province  ZZ)  show  a  smooth,  relatively  broad  distribution  of 
void  ratios. 

The  shelf  grain  densities  usually  lie  quite  close  to  a  value 

Q 

of  2.6  gm  cm”  while  slope  grain  densities  are  much  closer  to  2.7 

-3 

gm  cm  .  Province  BA  shows  a  large  percentage  of  grain  densities 
above  2.8  gm  cm-  --  an  unusually  high  value. 


t 


e.  Province  AC 


Province  BC 


Figure  2-7.  Grain  Density  Empirical  Cumulative  Distributions 
The  left-hand  curves  are  for  shelf  sediments, 
the  right-hand  curve  is  for  slope  sediments  of 
same  size  type. 
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Figure  2-7  (cont'd) 

Province  ZZ  represents  continental  terrace  (shell 
slope,  rise)  sediments  as  a  whole. 


e. 


Province  AC 


f.  Province  BC 


Figure  2-8. 


Mean  Specific  Surface  Empirical  Cumulative  Distri¬ 
bution  Functions.  The  left-hand  curves  are  for 
shelf  sediments,  the  right-hand  curve  is  for  slope 
sediments  of  the  same  grain  size  type. 


The  mean  specific  surface  shows  a  stronger  dependence  on 
sediment  grain  size  type  than  it  does  on  shelf/slope  differences. 
Sands  (AA  and  AB)  lie  in  the  range  0  to  10,000  cm--*-;  dirty  sands 
(AB  and  BB)  lie  approximately  between  5,000  to  15,000  cm-^ ;  sandy 
muds  (AC  and  BC)  in  the  range  5,000  to  30,000  cm“l;  silts  (AD  and 
BD)  in  the  range  10,000  to  20,000  cm-^;  and  muds  (AE  and  BE)  in 
the  range  20,000  to  60,000  cm-^.  Continental  terrace  sediments 
overall  (ZZ)  vary  smoothly  from  0  to  60,000  cm-*. 

Other  physical  property  characteristics  of  these  provinces 
are  not  available  from  the  PHYPR0SE  data  base,  yet  are  needed  as 
input  to  the  PHYSED  model.  Three  parameters  considered  important 
to  the  results^  are  frame  Poisson's  ratio,  frame  log  dec  of 
compress ional  waves,  and  frame  log  dec  of  shear  waves. 

We  define  ranges  for  these  three  properties  in  Table  2-9. 
The  values  tabulated  there  are  based  upon  guidelines  published  by 
Hamilton  .  Poisson's  ratio  for  marine  sediments  should  fall  in 
the  range  from  0.1  to  a  theoretical  upper  limit  of  0.5.  In 
sands,  values  between  0.1  and  0.3  are  typical,  while  in  muds, 
higher  values  are  expected.  Thus,  we  use  0.2  to  0.4.  The  frame 
log  dec  for  compressional  waves  is  expected  to  lie  in  the  range 
0.00  (for  negligible  loss  at  grain  to  grain  contacts)  to  0.3  in 
lossy  frames.  Silts  and  clays  are  expected  always  to  show  loss, 
so  their  range  is  assigned  as  0.1  to  0.2;  sands  are  allowed 
negligible  loss  and  assigned  a  range  of  0.0  to  0.15.  Shear  log 
decs  for  the  frame  are  expected  to  be  higher  than  compressional 
log  decs;  therefore,  ranges  are  assigned  as  in  Table  2-9.  It  is 
assumed  that  values  of  these  three  parameters  are  evenly  dis¬ 
tributed  across  the  ranges  given. 


Table  2-9.  Ranges  for  Poisson1 
Frame  Compress ional  Log  Dec  and  Frame 


2.5  Physical  Properties  Summary 

Based  on  our  survey  of  available  sediment  physical  property 
data,  we  find  texture  (grain  size)  information  ten  times  more 
readily  available  than  void  ratio/porosity  information  and  one 
hundred  times  more  readily  available  than  geoacoustic  infor¬ 
mation.  Based  on  the  data  from  233  sites  in  the  PHYPROSE  data 
base,  we  find  that  void  ratio/porosity  varies  not  just  with  sedi¬ 
ment  texture  but  also  with  depth  of  the  seafloor  and  ocean  area. 
We  therefore  define  provinces  based  on  grain  size  type  and  sea¬ 
floor  depth.  We  retain  the  option  of  introducing  ocean  area  as  a 
provincing  criterion  if  we  cannot  achieve  sufficiently  reduced 
variance  in  the  BL  predicted  for  a  province.  We  use  the  entire 
PHYPROSE  data  base,  with  its  data  at  4,355  different  sites,  to 
define  the  distributions  of  grain  density,  void  ratio/porosity, 
and  mean  specific  surface  in  the  provinces  on  the  shelf  and 
slope.  We  cite  the  literature  to  establish  the  ranges  of  three 
other  parameters  --  Poisson's  Ratio,  frame  compressional  log  dec, 
and  frame  shear  log  dec.  The  remaining  Biot/Stoll  physical 
inputs  are  presumed  fixed  and  therefore  independent  of  province. 


3.0  THE  GEOACODSTIC  PROFILES 


The  results  of  the  Biot/Stoll  calculations  are  the  geo¬ 
acoustic  profiles  of  speeds  and  attenuations  for  compress ional 
and  shear  waves  in  the  sediment.  From  these  profiles,  BL  will  be 
calculated  as  described  in  Section  4.  In  the  present  section, 
however,  we  describe  how  the  Biot/Stoll  computations  use  the  dis¬ 
tributions  generated  for  each  province,  and  we  show  the  distribu¬ 
tion  of  the  resulting  speed  and  attenuation  profiles.  We  also 
compare  the  results  of  the  computations  to  the  few  compressional 
speed  measurements  available  in  the  PHYPROSE  data  base. 

3.1  Monte  Carlo  PHYSED  Calculations 

The  many  options  available  for  defining  the  thirteen  PHYSED 
inputs,  as  illustrated  by  Figure  1-1  and  the  list  of  "genesis 
parameters"  in  Table  1-3,  were  reduced  to  those  few  which  could 
be  characterized  for  each  province.  PHYSED  software13  was 
modified  to  hardwire  the  input  generation  as  follows. 

The  fluid  properties  were  kept  fixed  for  all  provinces  at  a 

density  of  1.025  gm  cm-3,  a  bulk  modulus  of  2.384  x  1010  dynes 
-2 

cm  ,  and  a  viscosity  of  .018  poise.  (These  values  correspond  to 
a  fluid  compressional  speed  of  1525.1  m  s_1.)  The  grain  bulk 
modulus  was  fixed  at  4.2  x  10  dynes  cm-  ,  a  value  appropriate 
for  quartz  grains,  and  the  grain  densities  were  sampled  from  the 
PHYPROSE  distributions  as  described  later. 

The  structure  factor  was  fixed  at  1.25  for  all  provinces. 
Though  this  is  a  good  value  for  sands,®’32  a  value  of  3.0  has 
produced  good  companions  with  silts  and  clays.®’13  However,  to 


22 

Brunson,  B.A.,  1983,  "Shear  Wave  Attenuation  in  Unconsolidated 
Laboratory  Sediments,"  Ph.D.  Dissertation,  Oregon  State 
University. 


avoid  introducing  an  artificial  break  among  sediments  in  the  same 
province  (i.e.,  province  ZZ) ,  a  constant  fixed  value  of  structure 
factor  is  used. 


Porosity  of  the  frame  is  calculated  from  void  ratio,  sampled 
from  the  empirical  distribution,  using  equation  3-1. 


(3-1) 


where  B  is  porosity  (decimal) 
and  e  is  void  ratio. 


The  procedure  is  described  later  for  sampling  void  ratio  from  the 
distribution  based  on  PHYPROSE  data. 


Permeability  (k,  in  cm2)  and  pore  size  parameter  (a,  in  cm) 
are  both  derived  from  the  porosity  value  obtained  as  above  and 
the  mean  specific  surface  (SQ)  sampled  from  the  empirical 

.  O  O 

distributions.  According  to  the  Kozeny  Carman  relationship  u  for 
permeabi lity , 


5.0'S4 


(1-c) 


(3-2) 


The  pore  size,  following  Hovem  and  Ingram,  can  be  set 
equal  to  twice  the  hydraulic  radius  (ratio  of  volume  filled  with 
fluid  to  the  wetted  surface).  Hence, 
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Carman,  P.C. ,  1956,  Flow  of  Gases  Through  Porous  Media, 
Academic  Press,  New  York. 
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Hovem,  J.M.  and  G.D.  Ingram,  1979,  "Viscous  Attenuation  of 
Sound  in  Saturated  Sand,"  J.  Acoust.  Soc .  Amer.,  66,  pp  1807- 
1812. 
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The  frame  shear  modulus  (^b)  can  be  derived  as  a  function  of 
depth  from  the  porosity  and  grain  density  using  a  procedure 
described  by  Stoll2^  based  on  empirical  relationships  between 
stress  and  shear  modulus  determined  by  Richart,  et  a_l.2®  The 
vertical  stress,  (in  dynes  per  cm2),  at  a  depth,  Z  (in  cm),  in 
the  sediment  is  computed  by  integrating  the  buoyancy-reduced 
weight  of  the  overlying  material,  i.e., 


1 1 


o 


(Pr-of)g 


dz 


(3-4) 


where  g  is  gravitational  acceleration  (set  to  980  cm  s-2),  8  is 
porosity  (kept  constant  with  depth),  Pr  is  density  of  grain 
(constant  with  depth),  and  z  is  the  integration  variable 
representing  depth.  Given  t1}  the  average  stress,  (used  in 

the  Richart  equations),  is 


T 

O 


(3-5) 


where  T2  and  T3  are  the  horizontal  components  of  stress.  A  value 
for  tq  is  obtained  by  assuming  and  t3  each  is  equal  to  Tj  in 
clays,  and  that  each  is  half  of  Ti  in  sands.  The  Richart 

equation  for  sand  is 


1230  ( 2 . 97-e) 2  Js 
ub  "  (1+e)  *  lo 


(3-6) 


The  equation  for  clav  is 


1630  (2.9 7-e )  2  . 

“b  (1+e)  ° 


(3-7) 
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Stoll,  R.D.,  1977,  "Acoustic  Waves  in  Ocean  Sediments," 
Geophysics ,  42,  pp  715-725. 
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The  result  of  the  calculation  is  a  frame  shear  modulus  that 
increases  with  depth  even  though  porosity  (hence  e)  is  held  con¬ 
stant  . 


The  imaginary  part  (ub*,  dynes  cm-2)  of  the  frame  shear 
modulus  is  a  function  of  the  real  part  and  the  logarithmic 
decrement  (A^ f  dimensionless) 


* 


u 


(3-8) 


T 

Neither  frame  log  decs  nor  imaginary  parts  are  well  known  for  any 
of  the  provinces;  so,  a  range  is  defined  for  log  decs  according 
to  Table  2-9,  log  decs  are  chosen  randomly  from  that  range,  and 
is  calculated  using  Equation  3-8. 


The  frame  bulk  modulus  (Kb>  dynes  cm'2)  is  derived  from  the 
frame  shear  modulus  already  described  and  a  Poisson's  ratio  (Rp, 
dimensionless)  sampled  from  a  range  defined  for  each  province 
(see  Table  2-9).  The  equation  is 


K 


b 


2  (1+R  > 
_ E_  ub 

3  (1-2) Rp 


(3-9) 


The  derivation  of  the  imaginary  part  (Kb*)  of  the  frame  bulk 
modulus  is  based  on  the  real  part  and  the  logarithmic  decrement 
(Ar)  sampled  from  a  range  for  each  province  (Table  2-9).  The 
equation  in  effect  is  analogous  to  (3-8). 


K, 


K 


(3-10) 


As  mentioned  above,  six  properties  were  sampled  from 
distributions  defined  for  each  province;  the  six  are:  void 
ratio,  grain  density,  mean  specific  surface,  Poisson's  ratio, 
frame  compressional  log  dec,  and  frame  shear  log  dec.  Here  we 
describe  how  this  sampling  was  accomplished. 


The  empirical  cumulative  distribution  functions  for  each 
property  for  each  province  were  stored  as  computer  files  of 
property  value/probability  pairs.  Six  independent  sequences  of 
random  numbers,  distributed  evenly  between  0.0  and  1.0,  were 
produced  using  a  FORTRAN  subroutine  with  six  different  seeds. 
Each  random  number  can  be  converted  to  a  value  of  one  of  the 
properties  using  that  property's  empirical  cumulative  distribu¬ 
tion  curve,  and  associating  that  random  number  with  its  position 
on  the  probability  axis.  For  example,  if  the  random  number  is 
0.2,  then,  for  the  curve  depicted  in  Figure  2-6k,  a  probability 
of  0.2  is  associated  with  a  void  ratio  of  about  1.6.  The  six 
sequences  of  random  numbers  were  thus  used  to  generate  sequences 
of  the  six  properties  for  each  province,  using  a  computer  sub¬ 
routine  to  interpolate  between  the  probability/property  pairs  in 
the  computer  files.  The  resulting  sequences  of  property  values 
are  thus  distributed  as  if  sampled  from  populations  with  the  same 
statistics  as  depicted  in  Figures  2-6,  2-7,  2-8  and  Table  2-9. 
The  same  six  sequences  of  random  numbers  were  used  to  generate 
property  values  in  each  of  the  eleven  provinces. 

The  computer  software  MCPHYS,  a  FORTRAN  program  listed  in 
Appendix  A,  samples  the  six  properties  named  above  for  a  given 
province  and  computes  the  Biot/Stoll  solutions;  then,  repeats  the 
process  using  the  next  elements  of  the  six  sequences,  until  a 
specified  number  of  Biot/Stoll  solutions  (or  runs)  are  obtained. 
Using  MCPHYS,  we  made  fifty  runs  per  province  in  order  to  have 
sufficient  samples  to  simulate  a  reasonable  distribution  of 
PHYSED  results  associated  with  that  province. 

At  this  point  it  is  worth  noting  two  features  of  this 
procedure.  First,  the  six  properties  are  treated  as  if  they  are 
independent  of  each  other  because  they  are  derived  from  random 
number  sequences  with  different  seed  values.  Thus,  for  example, 
a  high  value  of  void  ratio  (resulting  from  a  random  number  near 
1.0)  can  occur  with  a  high,  low,  or  medium  value  of  any  of  the 
others  --  say  grain  density  --  because  a  different  random  number 
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sequence  is  used  with  that  other  property.  We  considered  the 
possibility  that  this  independence  might  not  be  realistic  in  some 
cases,  especially  with  respect  to  the  two  properties  void  ratio 
and  mean  specific  surface.  As  a  general  rule  sands  (larger 
grains  with  smaller  specific  surfaces)  tend  to  be  more  tightly 
packed  (hence  have  smaller  void  ratios).  So  one  might  suspect 
that  void  ratio  and  specific  surface  are  correlated  in  the 
sediments.  Using  the  data  in  PHYPROSE,  however,  we  were  unable 
to  demonstrate  a  strong  enough  correlation  (see  Figure  3-1  for 
Province  ZZ)  between  void  ratio  and  mean  specific  surface;  so  we 
considered  it  reasonable  to  treat  them  as  independent. 

The  second  feature  of  our  procedure  is  that  the  samples  of  a 
given  property  from  different  provinces,  but  the  same  position  in 
sequence,  are  correlated.  This  is  because  the  same  random  number 
sequence  is  used  for  a  given  property  for  all  provinces.  For 
example,  if  the  fourth  random  element  of  the  Poisson's  Ratio 
sequence  were  high  (say  .85)  for  province  AA,  than  it  would  be 
high  (  =  .85)  for  province  AB,  and  AC,  etc.  The  Poisson's  Ratio 
that  results  from  this  high  random  number  may  be  different  from 
one  province  to  the  next,  but  it  will  be  true  for  all  provinces 
that  the  fourth  Poisson's  Ratio  will  be  high  (in  the  limit  of 
many  samples,  higher  than  85%  of  the  Poisson's  Ratios  in  that 
province).  This  is  an  important  point  because  it  will  allow  us, 
when  comparing  provinces,  to  detect  small  differences  in  BL 
between  two  province.  We  can  do  this  by  accumulating  differences 
between  individual  members  of  the  two  provinces  that  occur  at  the 
same  position  in  sequence  (e.g.,  the  fourth  BL  of  province  AA 
subtracted  from  the  fourth  BL  of  province  AB) . 

3.2  Distribution  of  Geoacoustic  Profiles 

The  results  of  the  PHYSED  calculations  are  compressional 
and  shear  speeds  and  attenuations  as  functions  of  frequency.  In 
addition,  because  we  use  a  shear  modulus  that  increases  with 
pressure,  our  Biot/Stoll  results  vary  with  depth,  especially  in 
the  upper  10  meters.  We  performed  our  calculations  at  several 
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frequencies  and  depths  appropriate  to  ASW  applications,  and  list 
them  in  Table  3-1. 


Table  3-1.  Seven  Frequencies 

and  Six  Depths  at  which 

Compressional  Wave  and  Shear  Wave  Speed  and 

Attenuation  are  Calculated. 

Frequencies  (Hz) 

Depths  (m) 

100 

0 

200 

1 

400 

5 

800 

10 

1600 

50 

3500 

100 

10000 

The  values  of  a  geoacoustic  parameter  as  a  function  of  depth  or 
frequency  we  term  a  geoacoustic  profile.  As  a  result  of  the 
Monte  Carlo  runs,  we  have  produced  profiles  with  6  depths  for 
each  of  six  geoacoustic  parameters  (Table  1-2),  at  each  of  seven 
frequencies,  for  each  of  50  runs,  for  each  of  eleven  provinces. 
Thus  we  have  generated  138,600  values  to  describe  the 
distribution  of  the  geoacoustic  properties  in  these  provinces  -- 
too  many  to  list  in  this  report. 


We  have  devised  a  presentation  that  depicts  the  general 
character  of  our  results  more  simply.  We  display  (depth) 
profiles  of  compressional  speed  and  shear  speed  for  an 
intermediate  frequency  (since  these  are  only  weakly  frequency 
dependent),  and  we  display  compressional  attenuation  and  shear 
attenuation  as  functions  of  frequency  for  an  intermediate  depth 
(since  these  are  strongly  frequency  dependent).  Instead  of 
showing  all  fifty  depth  profiles  for  speeds  we  define  a  depth 
profile  of  mean  speeds,  and  a  square  root  of  variance  in  speeds 
at  a  given  depth.  In  addition  we  can  define  a  "typical"  profile, 
and  two  "extreme"  profiles.  These  six  profiles  then  describe  the 
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essentials  of  the  distribution  of  compressional  or  shear  speeds 
in  the  province.  In  a  similar  way,  six  frequency  profiles  can 
describe  the  essentials  of  the  distribution  of  attenuation  in  the 
province  for  either  compressional  or  shear  waves.  The  FORTRAN 
program  PROFIL,  listed  in  appendix  A,  computes  the  mean  and  vari¬ 
ance  of  any  profile  through  the  data  and  selects  the  typical  and 
extreme  profiles. 


In  Figures  3-2  through  3-5  we  present  some  of  our  results  in 
the  form  just  described.  Figure  3-2  shows  the  distribution  of 
compressional  speed  depth  profiles  in  three  provinces:  ZZ  (all 

continental  terrace  sediments  combined),  AA  (continental  shelf 
sands),  and  BE  (continental  slope  muds).  These  profiles  are  well 
behaved  in  the  sense  that  the  typical  profile  and  the  mean  are 
nearly  identical  in  location  and  shape,  and  even  the  extreme 

profiles  are  roughly  parallel  to  the  mean  profile.  Note  the 
differences  among  the  provinces  in  location  (compressional  speed 
of  the  mean)  and  spread  (compressional  speed  difference  between 
the  mean-plus-the-root-variance  and  the  mean-minus-the-root- 
variance)  of  the  profiles.  Shelf  sands  (AA)  have  the  highest 

speeds  at  about  1620  ms-^  while  slope  muds  have  the  lowest  at 
about  1480  ms-^.  The  provinces  are  well  separated  with  almost  no 
overlap  of  their  distribution  envelopes.  Province  ZZ  shows  an 

intermediate  speed  with  a  distribution  that  overlaps  both  of  the 
other  provinces.  It  is  precisely  these  kind  of  differences  that 
we  anticipate  will  lead  to  reduced  uncertainity  in  BL  predictions 
when  sediments  are  divided  into  provinces. 

In  Figure  3-4  we  illustrate  the  distribution  of  compression¬ 
al  attenuations  in  the  same  three  provinces:  ZZ,  AA,  and  BE. 

These  profiles  again  are  well  behaved;  but,  there  is  a  decided 
asymmetry.  Deviations  below  the  typical  profile  are  small  com¬ 
pared  to  deviations  above  the  it.  Again  note  the  differences 
among  the  provinces  in  location  and  spread  of  the  profiles.  Slope 
muds  have  lower  attenuations  and  smaller  variance  than  shelf 


Figure  3-2.  Compressional  speed  versus  depth  for  three 
provinces.  For  each  province  we  plot  the  mean  speed  at  a  given 
depth  (curve  1) ,  the  mean  speed  plus  the  square  root  of  the 
variance  (curve  2),  the  mean  speed  minus  the  square  root  of  the 
variance  (curve  3),  an  actual  profile  with  typical  values  (curve 
5),  and  two  actual  profiles  with  extreme  values  (curves  4  and  6). 
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Figure  3-5. 


Same  as  Figure  3-  4  but  for  Shear  Attenuation  versus 
Frequency  for  Three  Provinces. 
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sands;  continental  terrace  sediments  overall  are  intermediate  in 
mean  attenuation  and  attenuation  variance. 


In  Figures  3-3  and  3-5  we  present  similar  results  for  shear 
speeds  and  attenuations,  respectively,  for  the  same  three 
provinces.  We  will  not  launch  into  a  thorough  and  formal 
evaluation  of  the  differences  among  the  provinces  in  this  section 
for  one  reason.  The  geoacoustic  profiles  are  not  the  operation¬ 
ally  relevant  parameters  --  bottom  loss  (BL)  is.  Thus  we  will 
reserve  our  comparisons  until  after  we  have  described  the 
derivation  of  BL  in  Section  4  and  presented  those  distributions. 

3.3  Accuracy  of  Geoacoustic  Distributions 

We  assess  the  accuracy  of  the  distributions  calculated  for 
our  provinces  by  comparing  them  to  available  observations  of  geo¬ 
acoustic  parameters. 

The  PHYPROSE  data  base  contains  geoacoustic  measurements 
only  of  compressional  speed  at  various  depths  in  the  upper  two  or 
three  meters.  Because  these  are  the  only  measurements  of  either 
geoacoustic  parameters  or  BL  available  to  us,  they  represent  our 
only  opportunity  to  assess  the  accuracy  of  the  distributions  we 
have  attributed  to  our  provinces. 

We  are  comparing  measurements  made  under  shipboard  condi¬ 
tions  (room  temperature  and  atmospheric  pressure)  with  predic¬ 
tions  of  in  situ  conditions,  under  fixed  seawater  properties, 
produced  by  our  calculations.  To  eliminate  these  differences  in 
conditions  we  calculate  the  ratio  of  speed  in  the  sediment  to 
speed  in  seawater  at  the  temperature  and  pressure  of  the  measure¬ 
ment  or  calculation.  This  speed  ratio  is  then  insensitive  to  the 
variations  of  water  properties  and  reflects  only  the  properties 
of  the  sediment.  Our  calculations  are  based  on  seawater  with  a 
density  and  bulk  modulus  that  correspond  to  a  seawater  sound 
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speed  of  1525.1  ms-1;  so  all  our  compressional  speed  results  are 
divided  by  this  value.  The  PHYPROSE  data  base  archives  not  only 
speeds  measured  in  the  sediment  core  but  also  the  temperature  of 
the  measurement.  This  temperature  is  used  with  a  salinity  of  35 

parts  per  thousand  (ppt)  and  a  pressure  of  1  atmosphere  in  the 

97 

Wilson  sound  speed  equation  to  compute  the  seawater  compres¬ 
sional  speed  to  divide  into  that  measurement. 

Compressional  speed  measurements  are  available  for  four 
provinces.  In  Figure  3-6  through  3-9  we  plot  the  sound  speed 

ratio  profiles  calculated  for  each  province  together  with  the 
speed  measurements  in  that  province.  Except  for  one  suspicious 
measurement  in  province  AC,  all  observations  fall  within  two 

standard  deviations  of  the  mean  profile,  indicating  that  our 
definition  of  provinces  is  consistent  with  the  observations.  We 
do  not  attempt  to  compute  the  probability  that  the  mean  of  the 

observations  is  the  same  as  the  mean  of  the  predictions  because 
in  no  case  do  we  have  a  sufficient  number  of  independent 
observations.  Most  observations  of  compressional  speed  ratio  for 
a  province  are  collected  at  a  few  locations  in  close  proximity 
from  a  single  cruise.  They  do  not  span  the  range  of  sediments 
that  belong  to  a  province.  We  are  unable  therefore,  with  the 

present  data,  to  quantify  the  accuracy  of  our  distributions.  We 
can  say  that  a  comparison  with  data  does  not  show  any  problems 
with  the  distributions.  This  is  a  negative,  but  necessary, 
result . 
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4.0  BOTTOM  LOSS  DISTRIBUTIONS 


The  geoacoustic  profiles  of  Section  3  are  here  used  to 
calculate  the  Bottom  Loss  (BL).  The  resulting  distributions  of 
BL  in  each  province  are  presented  and  the  statistics  that  will 
serve  as  criteria  for  evaluating  improvement  are  tabulated. 

4.1  Calculation  of  Bottom  Loss 

In  the  bottom  loss  calculation  we  treat  the  sediment  as  a 
fluid  using  the  NORDA  computer  program  REFLEC^l.  While  most 
sediments  possess  some  rigidity  the  shear  wave  velocities  are 
small  enough  so  that  bottom  loss,  for  our  purposes,  will  not  be 
appreciably  influenced.  This  is  not  necessarily  true  for 
environments  more  complex  than  a  sediment  halfspace;  however,  the 
sediment  halfspace  (homogenous  sediment)  is  the  only  environment 
type  that  we  consider  in  this  study. 

Figure  4-1  shows  the  magnitude  of  the  reflection  coefficient 
|R|  at  100  Hz  for  a  clay  where  the  clay  is  modeled  as  a  fluid 
and,  again,  as  a  porous  viscoelastic  solid.  The  results  are 
essentially  identical.  Table  4-1  lists  the  geoacoustic  inputs  to 
the  model  for  the  clay  case.  These  inputs  were  smoothed  before 
running  the  reflection  coefficient  program. 

Table  4-2  lists  the  geoacoustic  inputs  to  the  models  for  a 
sand  case.  Figure  4-2  shows  IrI  at  100  Hz  plotted  for  sand  where 
the  sand  is  modeled  as  a  fluid  and,  again,  as  a  solid.  At  angles 
below  critical  (including  5°  where  we  choose  to  look  at  BL)  the 
solid  model  shows  higher  loss  due  to  excitation  and  subsequent 
attenuation  of  the  shear  waves.  The  geoacoustic  inputs  to  the 
model,  for  Table  4-2  which  were  smoothed,  show  high  shear 
velocities  to  provide  a  "worst-case"  input  set. 

When  we  model  sands  as  a  fluid,  then,  we  slightly 
underestimate  the  true  BL  value.  This,  however,  will  not 
appreciably  influence  our  analysis  of  BL  variance  because  the 


Table  4-1.  Geoacoustic  Model  for  Clay 


VP 

kP 

vs 

^s 

p  Depth 

(m/s) 

(dB/m/kHz) 

(m/s) 

(dB/m/kHz) 

(g/cm3)  (m) 

1470.69 

.000 

1470.91 

.000 

1471.08 

.000 

1471.38 

.000 

1472.21 

.000 

1473.54 

.000 

5.00 

11.13 

13.23 

15.74 

19.79 

23.53 


20.00 

15.60 

13.10 

11.03 


8.7' 

7.3 


L .  42 

0 

L .  42 

5 

L .  42 

10 

L .  42 

20 

L  .42 

50 

L .  42 

100 

Table  4-2.  Geoacoustic  Model  for  Sand 


VP 

kP 

Vs  ks  P 

D< 

(m/s) 

(dB/m/kHz) 

(m/s)  (dB/m/kHz)  (g/cm3) 

1636.69 

.3983 

87.76 

18.0 

2.02 

0 

1656.31 

.3764 

167.76 

14.48 

2.02 

5 

1664.45 

.3683 

199.51 

12.17 

2.02 

10 

1675.98 

.3577 

237.25 

10.24 

2.02 

20 

1698.88 

.3394 

298.33 

8.14 

2.02 

50 

1724.74 

.3226 

354.78 

6.85 

2.02 

100 

variance  is  principally  dependent  on  the  existence  of  the 
critical  angle,  which  is  in  turn  dependent  on  other  features  of 
the  sediment  than  conversion  to  shear  waves. 

4.2  Bottom  Loss  Distributions 

Figures  4-1  and  4-2  display  the  two  types  of  BL  curves 
expected  when  acoustic  waves  impinge  on  fluid  sediments  with 
densities  greater  than  seawater.  Figure  4-1  shows  the 
characteristic  curve  for  reflection  from  sediments  with 
compressional  speeds  lower  than  seawater.  There  is  an  angle  of 
intromission  (about  15°  in  Figure  4-1)  at  which  all  of  the 
incident  energy  is  transmitted  (hence  there  is  infinite  bottom 
loss  at  the  angle).  For  fixed  water  properties,  as  in  our  study, 
the  angle  of  intromission  will  vary  with  the  density  of  the 
sediment  and  the  compressional  speed  in  the  sediment.  At  grazing 
angles  from  the  angle  of  intromission  to  zero,  the  reflection 
coefficient  rises  steeply  to  one,  (and  BL  drops  steeply  to  zero) 
at  a  near-uniform  rate. 

Figure  4-2  shows  the  characteristic  curve  for  reflection 
from  sediments  with  compressional  speeds  higher  than  seawater. 
There  is  a  critical  angle  (about  32°  in  Figure  4-2)  below  which 
there  is  nearly  complete  internal  reflection  (and  near-zero  BL). 
Because  of  attenuation  in  the  sediment  actual  BL  at  these  angles 
can  vary  from  0  dB  to  3  dB.  At  the  critical  angle,  reflection 
drops  drastically  (and  BL  rises  rapidly);  and,  at  grazing  angles 
a  few  degrees  greater  than  the  critical  angle,  BL  levels  out. 
The  critical  angle  is  a  function  of  the  compressional  speed  in 
the  sediment. 

The  above  characteristics  of  BL  curves  for  sediments  will  be 
important  in  interpreting  our  results.  The  BL  at  a  grazing  angle 
of  5°  is  considered  indicative  of  the  lower  angles  important  to 
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longer  range  propagation  in  ASW  applications.  Thus  we  calculate 
and  present  distributions  of  BL  at  5°  to  characterize  our 
provinces.  Note  that  5°  is  typically  in  the  near-zero  BL  range 
for  fast  sediments  (i.e.,  sediments  in  which  compressional  sound 
speed  is  greater  than  in  seawater);  and  that  5°  is  typically  in 
the  steep  BL  ramp  between  0°  and  the  angle  of  intromission  for 
slow  sediments  (in  which  compressional  speed  is  less  than  in 
seawater) . 

4.3  BL  in  Provinces  AA  and  BE  and  ZZ 

Figure  4-3  shows  histograms  of  BL  at  100  Hz  for  three 
provinces:  AA  (shelf  sands),  BE  (slope  muds),  and  ZZ  (all 
continental  terrace  sediments  combined).  Figure  4-4  shows 
histograms  of  BL  for  the  same  three  provinces  but  at  1600  Hz. 

Appendix  B  contains  twenty-two  BL  histograms  that  show  the 
same  two  frequencies  in  each  of  the  eleven  provinces  of  this 
study.  All  these  histograms  have  a  BL  resolution  of  0.5  dB. 

Figure  4-3a,  for  province  AA,  shows  a  BL  distribution 
expected  for  fast  sediments.  All  values  lie  between  0  and  about 
3  dB.  This  is  consistent  with  the  compressional  speeds  presented 
for  the  province  in  Figure  3-2  because  they  are  everywhere 
greater  than  the  seawater  speed  of  1525  m  s-1.  In  spite  of  the 
wide  range  of  speeds  associated  with  the  province,  the  BL  values 
are  distributed  tightly  near  0  dB.  There  is  a  BL  mode  in  the  0 
to  0.5  dB  interval,  the  mean  BL  is  .8  dB,  and  the  root  variance 
is  less  than  1  dB. 

Figure  4-3b,  for  province  BE,  shows  a  very  different 
distribution.  The  BL  values  are  higher  and  are  distributed  more 
widely  with  a  mode  in  the  7.5  to  8.0  dB  interval,  a  mean  of  8.8 
dB,  and  a  root  variance  of  almost  3  dB.  This  is  consistent  with 
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Figure  4-4.  Histograms  (with  0.5  dB  resolution)  representing  the 
distribution  of  Bottom  Loss  at  1600  Hz  in  three  provinces:  shelf 
sands  (AA) ,  slope  muds  (BE) ,  and  all  continental  terrace  sediments 
combined  (ZZ) .  Total  number  of  BL  values  in  each  province  is  50. 


slow  sediments  and  a  angle  of  intromission  varying  about  a  value 
like  that  shown  in  Figure  4-1.  Province  BE  is  definitely  a  slow 
sediment  because,  as  shown  for  Figure  3-2,  its  compressional  wave 
speeds  are  always  less  than  the  seawater  value  of  1525  m  s-1. 

Province  ZZ,  whose  BL  histogram  is  Figure  4-3c,  shows  two 
modes;  one  near  0  dB  associated  with  the  fast  sediments  on  the 
terrace  and  one  near  8  dB,  apparently  associated  with  the  slow 
sediments.  The  mean  BL  for  province  ZZ  (6.2  dB)  falls  between 
the  means  for  province  AA  and  BE.  The  root  variance  is  greater, 
at  3.3  dB,  than  in  either  province. 

The  story  is  essentially  the  same  at  1600  Hz,  as  told  in 
Figure  4-4.  These  three  provinces  clearly  demonstrate  an 
advantage  to  provincing.  By  the  criteria  set  up  in  Section  1.2 
(page  11),  we  have  shown  that: 

•  Provinces  defined  on  the  basis  of  water  depth  and  grain 
size  are  based  on  information  about  a  hundred  times  more 
plentiful  than  geoacoustic  and  BL  measurements  (as  shown 
in  Section  2.1) . 

.  .  .and. . . 

•  The  variances  found  for  province  AA  and  BE  are  lower  than 
those  for  ZZ. 

•  The  means  found  for  provinces  AA  and  BE  are  different 
those  for  ZZ. 


This  is  our  definition  of  improved  capability;  the  provinces  AA 
and  BE  satisfy  this  definition  qualitatively.  In  Section  5.0  we 
present  statistical  tests  to  determine  whether  these  differences 
in  means  and  variances  are  significant. 


4.4  BL  in  other  Provinces 

The  distributions  of  BL  in  the  other  provinces  (appendix  B) 
are  not  always  as  simple  as  those  already  presented.  There  is 
some  overlapping  of  results,  and  some  high  variances.  Problems 
arise  that  can  be  associated  with  deficiencies  in  the  contents  of 
the  PHYPROSE  data  base.  However,  we  believe  the  data  base 
problems  can  be  alleviated  by  the  addition  of  more  representative 
data,  and  that  the  overlapping  of  results  can  be  reduced  by 
redefining  provinces  --  both  of  which  are  beyond  the  scope  of 
this  study.  As  we  show  in  the  next  section,  significant 
improvement  is  attained  by  the  provincing  presented  here; 
consequently,  there  is  sufficient  justification  within  this  study 
to  pursue  this  approach  in  future  efforts. 

Provinces  AB  and  AC  (shelf  mixtures  with  sands  --  dirty 
sands  and  sandy  muds)  both  show  results  analogous  to  Province  AA 
at  100  Hz.  They  have  strong  modes  at  0  dB  BL  and  root  variances 
less  than  1  dB,  indicative  of  fast  sediments.  However,  at  1600 
Hz  both  provinces  show  a  few  cases  (about  10%)  with  exceedingly 
large  BL  values  (above  20  dB)  .  These  deviates  displace  the 
distribution  mean  away  from  the  mode,  and  introduce  two  of  the 
largest  root  variances  encountered  in  our  study  (8.3  dB  and  6.3 
dB).  These  root  variances  are  double  the  root  variance  of 
continental  terrace  sediments  as  a  whole  (Province  ZZ  with  3.9 
dB) . 

In  each  of  these  provinces  at  1600  Hz,  the  compressional 
speed  distribution  has  a  mean  value  about  one  standard  deviation 
above  the  seawater  speed  (see  Figure  3-7).  Thus,  in  about  15%  of 
the  cases,  AB  and  AC  should  behave  as  slow  sediments  and  give  BL 
values  in  excess  of  3  dB.  These  slow  speeds  are  found  only  in 
the  upper  10  meters  of  provinces  AB  and  AC;  consequently  waves  at 
100  Hz,  responding  to  average  conditions  in  the  upper  fifty  to 
one  hundred  meters,  do  not  behave  as  if  interacting  with  a  slow 
sediment . 


vv 


There  still  remains  the  question  of  why  BL  associated  with 
the  slow  sediment  in  the  provinces  is  so  high.  The  maximum  BL 
encountered  in  the  slow  sediments  of  province  BE  is  15  dB;  yet 
provinces  AB  and  AC  show  BL  values  over  20  dB.  This  difference 
is  reasonable  given  the  greater  density  of  the  sandy  sediments  of 
AB  and  AC  over  the  clay  sediments  of  province  BE.  The  higher 
density  sediments  have  angles  of  intromission  at  lower  grazing 
angles  than  the  lower  density  sediments;  hence,  a  grazing  angle 
of  5°  is  closer  to  total  loss  in  sandy  sediments  than  in  clayey 
sediments . 

The  following  argument  makes  this  case  more  rigorously.  The 

definition  of  the  angle  of  intromission  (©T,  relative  to  normal 

28  ^ 

incidence)  is  given  by  Kinsler,  e_t  al.  in  their  equation  6.34 


sm  6  =  W - T 

1  V  1- (d1/p2) 


(4.1) 


where:  the  subscript  1  refers  to  seawater  and  the  subscript  2 

refers  to  the  sediment 
-  is  density 

r  is  impedance,  equal  to  the  product  of  compressional 
speed  and  density 

and  I  is  measured  in  degrees  from  the  perpendicular  to  the 

sediment  surface 


In  our  notation,  this  becomes 


cos  0  = 

I 


l-(Pf  V£/o  V  ) 


1“ (pf  /  0] 


(4.2) 


28fCinsler,  L.E.,  A.R.  Frey,  A.B.  Coppens,  and  J.V.  Sanders,  1980, 
Fundamentals  of  Acoustics,  3rd  Ed.,  John  Wiley  &  Sons  Inc.,  New 
York,  New  York. 
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where  Pf=  density  of  fluid  =  1.025  gm  cm 

Vf  =  compressional  speed  in  fluid  =  1525  m  s-1 

P  =  sediment  density  —  a  function  of  grain  density  and 
porosity 

Vp  =  compressional  speed  in  sediment 

0  =  grazing  angle  of  intromission  (degrees  from  the 

1  sediment  surface) 

We  use  equation  4.2  to  show  that  the  angle  of  intromission 

0T  is  much  smaller  for  province  AB  than  for  BE.  Typical  values 
in  AB  under  slow  conditions  for  speed  and  density  are  1510  m  s 
and  1.86  gm  cm-3  (based  on  a  void  ratio  of  1.0  and  a  grain 
density  of  2.70).  In  province  BE  all  conditions  are  slow  and  a 
typical  speed  is  1480  m  s~^,  and  a  typical  density  is  1.43  gm  cm-3 
(based  on  a  void  ratio  of  3.0  and  a  grain  density  of  2.65). 
These  values  lead  to  a  grazing  angle  of  intromission  of  5.4°  for 
AB  and  14.8°  for  BE.  It  is  clear  that  our  angle  of  5°  is  near 
total  loss  in  the  sandy  provinces.  Thus  in  those  few  cases  of 
slow  sediments  in  provinces  AB  and  AC,  the  higher  densities  lead 
to  smaller  angles  of  intromission  and  therefore  higher  BL  values 
than  are  found  in  province  BE. 

The  pattern  introduced  by  the  distribution  of  BL  at  1600  Hz 
in  provinces  AB  and  AC  is  repeated  and  reinforced  through  most  of 
the  remaining  provinces.  Two  processes  control  the  bottom  loss 
within  a  province,  one  at  sediment  compressional  speeds  greater 
than  those  of  seawater  (characterized  by  low  bottom  loss)  and  one 
at  speeds  lower  than  those  of  seawater  (characterized  by  high 
bottom  loss).  The  competition  of  these  two  processes  leads  to 
distributions  that  are  bimodal  for  many  of  the  provinces.  We  can 
define  a  low  BL  mode  as  that  BL  bin  below  4  dB  with  more  counts 
than  any  other  bin.  We  can  define  a  high  BL  mode  as  that  BL  bin 
above  4  dB  with  more  counts  than  any  other  bin.  We  select  no 
mode  unless  there  are  at  least  two  counts  in  a  bin;  and  we  name  a 
bin  by  the  BL  value  at  the  lower  end  of  its  interval.  With  these 
guidelines  <e  are  able  to  construct  Table  4-3,  which  shows  the 
modes  in  the  BL  distributions  for  100  Hz,  and  Table  4-4  for  1600Hz 


Table  4-3.  Modes  in  the  BL  Distributions  for  100  Hz 


sand 

dirty 

sands 

sandy 

muds 

silt 

mud 


Shelf 

Slope 

province 

name 

low  BL 
mode 

high  BL 
mode 

low  BL 
mode 

high  BL 
mode 

province 

name 

AA 

0.0 

- 

0.0 

15.5 

BA 

AB 

0.0 

- 

0.0 

11.0 

BB 

AC 

0.0 

- 

0.5 

9.0 

BC 

AD 

0.0 

11.5 

2.0 

12.0 

BD 

AE 

0.0 

9.5 

0.0 

7.5 

BE 

Table  4-4.  Modes  in  the  BL  Distributions  for  1600  Hz 


N.  water 

Ny  depth 

grain's 
size  N. 

class  \ 

Shelf 

Slope 

province 

name 

low  BL 
mode 

high  BL 
mode 

low  BL 
mode 

high  BL 
mode 

province 

name 

sand 

AA 

0.0 

- 

0.5 

BA 

dirty 

sands 

AB 

0.0 

- 

0.0 

11.5 

BB 

sandy 

mud 

AC 

0.0 

15.0 

0.5 

9.0 

BC 

silt 

AD 

0.0 

10.5 

- 

13.5 

BD 

mud 

AE 

0.0 

9.5 

- 

7.0 

BE 

These  tables  show  a  tendency  for  the  high  mode  to  occur  at  BL 

values  that  decrease  from  shelf  to  slope  and,  especially,  from 
sand  to  mud.  This  trend  is  driven  by  the  effects  of  void  ratio 
on  the  angle  of  intromission  via  sediment  density. 

Of  course,  not  all  modes  are  equal,  and  in  several  provinces 
one  mode  dominates.  Where  two  modes  compete  the  variance  is 
high;  where  one  mode  dominates  the  variance  is  low.  Only 

province  AB,  at  1600  Hz,  shows  high  variance  without  being  able 

to  define  a  high  BL  mode.  Using  a  root  variance  of  4  dB  to 

separate  2-mode  provinces  with  a  dominant  mode  from  2-mode 
provinces  with  competing  modes,  we  can  make  Table  4-5.  The 
provinces  with  competing  modes  do  not  show  improvement  over  a 
province  that  groups  all  continental  terrace  sediments  together. 

Table  4-5. 


Province 

100  Hz 

1600  Hz 

AA 

AB 

AC 

competing  modes 

AD 

competing  modes 

competing  modes 

AE 

competing  modes 

competing  modes 

BA 

competing  modes 

BB 

competing  modes 

competing  modes 

BC 

competing  modes 

competing  modes 

BD 

BE 

The  sand  sediments  on  the  slope,  province  BA,  do  not  give 
a  BL  distribution  as  similar  to  province  AA  as  we  expected. 
There  are  many  observations  of  high  BL  (about  50%  of  the 
observations  have  BL  greater  than  3  dB)  indicating  a  significant 
fraction  of  slow  sediments  in  the  province.  This  is  not 
considered  reasonable.  Upon  close  inspection  we  found  that  the 


speeds  below  seawater  compressional  speeds.  Thus  the  anomalous 
grain  densities  created  a  spurious  high  BL  mode.  Province  BA  is 
unreliable  because  the  data  in  PHYPROSE  for  this  province 
violated  our  assumption  that  they  were  representative.  This  kind 
of  error  can  only  be  alleviated  by  increasing  the  amount  of  data 
in  PHYPROSE  before  pursuing  further  evaluations  of  provincing. 
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5.0  RESULTS 

5.1  Statistical  Comparison  of  BL  Among  Provinces 

The  trends  and  distinctions  discussed  in  Section  4  are 
informative  but  are  not  associated  with  a  level  of  significance 
or  a  degree  of  confidence.  For  these  measures  of  reliability  we 
turn  to  analysis  of  variance  (ANOVA)  procedures.29 

Table  5-1  lists  the  quantities  used  in  an  analysis  of 
variance  for  BL  in  all  eleven  provinces  at  100  Hz,  and  Table  5-2 
does  the  same  at  1600  Hz. 

At  100  Hz  the  F  statistic  is  23.1  and  at  1600  Hz  the  F 
statistic  is  31.1,  which  indicate  highly  significant  differences 
in  the  means  associated  with  the  provinces  at  both  frequences. 
An  F  statistic  of  greater  than  2.7  is  significant  at  the  0.5% 
level.  Even  though  the  means  have  no  particular  physical  meaning 
in  the  provinces  where  the  distributions  are  bimodal,  and  the 
analysis  of  variance  depends  on  near  normal  distributions,  the 
extremely  high  F  statistics  do  indicate  significant  differences 
in  the  distributions  as  reflected  by  the  means. 

We  employ  the  Q  statistic29  to  isolate  those  provinces  whose 
means  are  different  from  each  other  at  the  95%  confidence  level. 
For  11  means  and  539  degrees  of  freedom,  Qq.05  ^Siven  by  Snedecor 
and  Cochran's  Table  A  1529)  is  less  than  4.70.  This  value  is 
multiplied  by  the  standard  error  of  the  mean,  s-  (the  square  root 
of  the  mean  square  within  the  province,  s  ,  divided  by  the  number 
of  samples  in  the  province,  n)  to  given  the  significant 
difference,  D,  between  means.  That  is 


D  =  Q 


0.05 


Sx  Q0.05 


2®Snedecor,  G.W. ,  and  W.G.  Cochran,  1973,  Statistical  Methods, 
Sixth  Edition,  Iowa  State  University  Press,  Ames,  Iowa. 
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2.69  dB. 


At  100  Hz  s2  =  16.3927  (Table  5.1)  so  D  = 


This 


indicates  that  the  means  of  provinces  AA,  AB,  and  AC  are  each 
significantly  different  from  (i.e.  less  than)  all  the  other 
provinces  (but  not  from  each  other).  The  only  other  significant 
difference  by  this  test  is  that  the  mean  of  province  AD  is  less 
than  the  mean  of  province  BB.  At  1600  Hz,  s2  is  31.80  so  D  is 
3.75.  This  leads  to  more  complicated  pairings  of  significant 
differences;  but  because  we  are  most  concerned  with  differences 
from  province  ZZ  (all  continental  terrace  sediments),  we  can 
reduce  the  complication.  Only  provinces  AA  and  AC  have  means 
significantly  less  than  the  mean  of  ZZ.  Only  provinces  BB,  BC, 
and  BD  have  means  significantly  greater  than  the  mean  of  AA. 


5.2  Separate  Analysis  of  High  and  Low  BL  Regimes 

Even  the  blind  analysis  of  variance  applied  in  Section  5.1 
is  able  to  ascertain  the  difference  between  provinces  with  a 
single  low  BL  mode  (AA,  AB,  and  AC  at  100  Hz)  and  provinces  with 
a  second  mode  at  high  BL  (most  of  the  remaining  provinces).  In 
this  section  we  wish  to  exploit  our  understanding  of  the  two 
regimes  (our  "model"  of  the  process)  to  make  more  subtle 
statistical  distinctions  --  such  as  the  difference  between  the 
high  BL  values  of  different  provinces.  The  low  BL  regime  occurs 
whenever  the  effective  compressional  speed  of  the  sediment  is 
greater  than  the  compressional  speed  of  seawater.  BL  values  are 
almost  always  less  than  3  dB  in  this  regime.  The  high  BL  regime 
occurs  whenever  the  effective  sediment  compressional  speed  is 
less  than  that  of  seawater.  In  this  regime  BL  values  are  almost 
always  greater  than  4.  We  will  now  generate  two  new 
distributions  from  each  histogram  in  Appendix  B;  one  low  BL 
distribution  from  0  to  3.5  dB,  and  one  high  BL  distribution  from 
3.5  to  40  dB.  We  will  do  comparisons  between  provinces  only  in 
similar  regimes;  e.g.,  compare  the  high  BL  distribution  of 
province  AA  with  the  high  BL  distribution  of  province  AB. 


Table  5-3  give  quantities  of  ANOVA  for  low  BL  conditions  at 
100  Hz,  while  Table  5-4  gives  those  quantities  for  high  BL 
conditions  at  100  Hz.  The  low  BL  (Table  5-3)  means  vary  from  .21 
dB  in  province  AB  to  2.23  dB  in  province  BD,  with  an  F  statistic 
of  17.17.  For  324  observations  spread  over  eleven  classes,  an  F 
statistic  over  2.71  is  significant  at  better  than  the  0.5%  level; 
therefore  there  are  highly  significant  differences  in  the  means. 
Again  employing  the  Q  statistic,  we  define  a  significant 
difference,  D.  The  equation  is  modified  from  Equation  5.1 
because  the  number  of  samples  in  a  province,  n,  is  now  a  function 
of  province;  according  to  reference  29,  page  278  the  standard 
error  of  the  mean  is  now  such  that 


D  =  Q 


0.05 


/f2(i  +  I 

2  \ni  nk 


(  5.2 


where  n^  is  the  number  of  samples  in  province  i  and 

nk  is  the  number  of  samples  in  province  k 

For  provinces  with  50  samples  D  =  .53  dB.  For  provinces  with  at 
least  18  observations  to  be  different  from  ZZ  (with  12 
observations),  D  must  be  greater  than  .78  dB.  There  are  many 
different  pairings  that  can  be  tested,  and  these  will  be 
summarized  and  plotted  later.  Here  suffice  it  to  say  that  at  100 

Hz  the  low  BL  means  in  provinces  AB  and  AC  are  significantly  less 

than  in  province  ZZ,  and  BD  is  significantly  greater  than  ZZ, 
with  95%  confidence. 

One  can  also  test  whether  there  are  significant  differences 
among  the  variances  for  the  provinces,  using  the  M/C  statistic 
(reference  29,  page  296).  This  value  was  computed  for  the  low  BL 
data  at  100  Hz  and  is  202.9  with  10  degrees  of  freedom.  Such  a 
result  indicates  real  differences  in  variance  at  better  than  the 
0.5%  level  of  significance.  Several  provinces  have  significantly 
lower  variance  than  ZZ,  as  will  be  shown  later,  and  no  province 
has  a  significantly  higher  variance. 


Table  5-4.  Analysis  of  Variance  Among  Provinces  for  Low  BL  at  5°  Grazing  Angle  at  1600  Hz 


Table  5-4  shows  similar  results  for  1600  Hz  low  BL 
conditions,  except  that  three  provinces  (BC,  BD,  BE)  had 
insufficient  data  (less  than  5  observations)  at  low  BL  to  be 
included  in  our  analysis.  Among  the  eight  remaining  provinces 
there  exists  an  F  statistic  of  10.19  which  is  significant  at 
better  than  the  0.5%  level  (greater  than  3.09  when  there  are  204 
observations  over  8  provinces).  The  M/C  statistic  is  47.74  with 
7  degrees  of  freedom  indicating  real  differences  in  the  variances 
among  the  provinces  at  the  0.5%  level  of  significance. 

Tables  5-5  and  5-6  show  the  results  of  ANOVA  calculations 
for  the  high  BL  cases  at  100  Hz  and  1600  Hz,  respectively.  Under 
high  BL  conditions  at  100  Hz  the  three  provinces  AA,  AB  and  AC 
have  insufficient  data  for  analysis  (less  than  5  observations). 
The  remaining  eight  give  an  F  statistic  of  6.77  which  is 
significant  at  the  0.5%  level  also.  (The  0.5%  level  is  defined 
to  be  under  3.09  for  218  observations  over  8  provinces).  The  M/C 
statistic  is  468.6  with  7  degrees  of  freedom  —  significant  at 
well  beyond  the  0.5%  level  of  chi  squared,  20.28.  Highly 
significant  differences  in  means  and  variances  therefore  occur 
among  the  provinces  at  100  Hz  in  the  high  BL  class. 

Finally,  at  1600  Hz  under  high  BL  conditions,  ANOVA  results 
are  summarized  in  Table  5-6.  Province  AA  has  too  few  observation 
to  be  included.  The  F  statistic,  at  11.53,  is  well  beyond  the 
0.5%  level  of  significance  value  of  2.60  based  on  320 
observations  over  10  provinces.  The  M/C  statistic,  at  410.8,  is 
well  beyond  the  0.5%  level  of  significance  for  chi  squared  with  9 
degrees  of  freedom,  23.59.  Again,  highlv  significant  differences 
in  means  and  variance  occur  among  the  provinces. 

Under  all  conditions,  100  Hz  and  1600  Hz,  low  BL  and  high 
BL,  the  provinces  defined  in  this  study  lead  to  significant 
differences  in  means  and  variances.  We  now  wish  to  indicate 
which  provinces  produce  those  differences.  The  results  of  our 
study  are  summarized  in  Table  5-7  for  100  Hz  and  Table  5-8  for 
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1600  Hz.  For  each  province  we  list  the  probability  that  BL  will 
be  less  than  3.5  dB.  If  sufficient  data  are  available  (5 
observations  or  more),  we  then  provide  the  mean  BL  and  variance 
for  the  low  BL  conditions,  and  the  mean  BL  and  variance  for  the 
high  BL  conditions.  We  also  provide  error  estimates  for  each  of 
the  values  in  Tables  5-7  and  5-8  to  aid  in  selecting  differences 
that  are  likely  to  be  significant. 

The  error  estimates  are  made  as  follows.  For  the  proportion 
of  observations  less  than  3.5  dB,  the  binomial  distribution  gives 
the  error.  The  99%  confidence  interval  is  encompassed  by  an 
error  in  sample  properties,  ££  given  by: 


where  p  is  the  proportion  of  the  sample  less  than  3.5  dB  and  n  is 
the  number  of  observations  in  the  sample  --  in  our  case,  always 
50. 


interval  based  on 

E  in  the  sample 
x 

5.4 

where  x  is  the  sample  mean 

s  is  the  sample  root  variance 
and  n  is  the  number  of  samples. 


For  the  mean  we  use  the  99%  confidence 

Student's  t  distribution  to  define  an  error 
mean. 

O 

£ 


E-  =  t 

x  .  01 


For  the  variance,  the  error  is  sensitive  to  the  kurtosis  of 
the  parent  distribution  which  we  estimate  from  the  sample 

distribution.  We  then  use  the  expected  error,  So  of  the  sample 

o  ®  29 

variance,  s  ,  from  a  distribution  with  kurtosis  Y. 


94 


(5.5) 


S„2=  a' 


|i{  1  +  ill  5-}| 


o  2 

where  az=  population  variance,  approximated  by  s 

f  =  population  kurtosis,  approximated  by  the  sample 

kurtosis 

f  =  degrees  of  freedom  =  n-1 


We  define  a  province  root  variance  upper  bound,  sU)  and  lower 
bound,  sl  as  follows. 


52  +  «  2 


s  .  es2 


(5.6) 


Finally,  we  define  the  error  estimate  of  the  root  variance 
£g  as  the  larger  of  the  differences  between  s  and  sUt  and  s  and 


Sg  =  maximum  (s^-s,  s-s^) 


The  results  of  Table  5-7,  for  100  Hz,  are  displayed  in 
Figures  5-1  through  5-5.  Figure  5-1  shows  the  proportion  of  low 
BL  conditions.  Only  provinces  AA,  AB,  AC  and  AD  have 
significantly  greater  proportions  (hence  probabilities)  of  low  BL 
conditions  than  occur  in  ZZ.  Additionally,  in  provinces  BA  and  BD 
it  is  likely  that  proportions  of  low  BL  conditions  are  greater 
than  in  ZZ,  but  the  confidence  level  of  this  distinction  is  not  as 
great  as  it  is  for  the  first  four  provinces. 


rs  represent  99%  confidence 


100  Hz 


Figure  5-2  shows  the  mean  BL  under  low  BL  conditions  at  100 
Hz.  While  several  significant  distinctions  can  be  made  between 
pairs  of  provinces,  only  province  AB  has  a  significantly  lower 
mean  BL  than  ZZ,  and  only  province  BD  has  a  significantly  higher 
mean  BL  than  province  ZZ.  A  good  part  of  the  reason  for  having 
only  two  provinces  distinct  from  ZZ  is  the  large  uncertainty  in 
the  ZZ  mean  because  of  the  small  sample  size  below  3.5  dB  in  that 
province. 

Figure  5-3  shows  the  variance  in  BL  under  low  BL  conditions. 
Provinces  AB,  AC  and  BD  have  significantly  lower  variances  than 
ZZ;  and,  province  BB  has  a  significantly  higher  variance.  The 
occurrence  of  significantly  lower  variances  in  some  provinces  is 
a  critical  factor  in  our  conclusion  that  there  is  an  advantage  to 
provincing  on  grain  size  and  water  depth. 

Figure  5-4  shows  the  mean  BL  at  100  Hz,  under  high  BL 
conditions.  Note  that  provinces  AA,  AB  and  AC  had  too  few 
observations  to  calculate  a  mean.  In  a  real  sense,  as  reflected 
in  Table  5-1,  these  three  provinces  are  significantly  different 
from  AA.  Among  the  remaining  provinces,  AE  and  BB  have 
significantly  higher  means  than  does  ZZ.  No  province  has  a 
significantly  lower  mean.  In  terms  of  the  mean  of  high  BL 
conditions  at  100  Hz,  province  BE  is  nearly  identical  in  location 
and  spread  to  province  ZZ.  This  property  was  first  mentioned  in 
Section  4.2.1  when  discussing  Figure  4-3.  The  significantly 
different  mean  values  in  some  provinces  is  another  of  the 
critical  factors  leading  to  our  conclusion  that  there  is  an 
advantage  to  this  provincing. 

Figure  5-5  shows  the  variance  of  BL  at  100  Hz  under  high  BL 
conditions.  Provinces  AA,  AB  and  AC  again  are  different  by 
virtue  of  having  no  high  BL  conditions.  Among  the  remaining 
provinces,  all  but  AE  and  BE  have  significantly  greater  variances 
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than  ZZ.  By  this  measure,  provincing  offers  no  advantage.  Once 
again  province  BE  is  most  like  ZZ. 

The  results  at  1600  Hz  are  qualitatively  similar  to  the 
findings  at  100  Hz.  The  1600  Hz  results  are  listed  in  Table  5-8 
and  plotted  in  Figures  5-6  to  5-10.  Provinces  AA,  AB,  AC  and  BA 
have  significantly  greater  proportions  of  low  BL  conditions  than 
does  ZZ  (Figure  5-6).  There  are  significant  differences  in  low 
BL  means  among  the  provinces;  but,  because  of  a  wide  uncertainty 
in  province  ZZ,  no  province  is  significantly  different  from  ZZ 
(Figure  5-7).  Provinces  AA,  AB,  AE  and  BB  have  significantly 
lower  variances,  under  low  BL  conditions,  than  does  ZZ,  and  no 
province  has  higher  variance  (Figure  5-8).  Under  high  BL 
conditions,  provinces  AC,  AD,  BB,  BC  and  BD  have  significantly 
higher  mean  BL  than  does  ZZ;  and  province  BE  is  nearly  identical 
to  province  ZZ  in  location  and  spread  (Figure  5-9).  Finally, 
several  provinces  have  significantly  larger  variances,  under  high 
BL  conditions,  than  does  ZZ  (Figure  5-10). 

5.3  Evaluation  of  PHYSED  Provincing 

We  restate  our  objective  (Section  1.2)  using  terminology 
introduced  in  this  report.  Use  of  the  PHYSED  model  can  be 
considered  an  improvement  in  capability  to  predict  BL  on  the 
continental  terrace  if  it  could  province  by  utilizing  information 
more  readily  available  than  geoacoustic  and  BL  measurements 
combined,  and  if  at  least  one  of  the  following  occurs. 

1)  The  variance  in  a  province  is  significantly  less  than 
the  variance  of  continental  terrace  sediments  combined 
( province  ZZ) . 

2)  The  mean  BL  in  a  province  is  significantly  different 
than  the  mean  of  province  ZZ. 
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1600  Hz 


Figure  5-6.  Significant  differences  between  ZZ  and  the  other  provinces  wrt  proportion 
of  low  BL  cases,  1600  Hz.  Bars  represent  99%  confidence  limits. 
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3)  The  probability  that  a  measurement  of  compressional 
speed  is  not  from  the  province  associated  with  it 
should  be  less  than  one  percent. 

The  provinces  considered  in  this  study,  useful  for  the 
PHYSED  model,  are  based  on  readily  available  information  for  a 
given  site  --  water  depth  and  qualitative  size  class.  Water 
depth  is  well  enough  known  worldwide  to  assign  a  location  to  the 
shelf  (0-200  m) ,  slope  (200  -  2000  m) ,  or  rise  (2000  -  4000  m) 
provinces.  As  shown  in  Section  2,  qualitative  size  class 

(texture)  information  is  on  the  order  of  one  hundred  times  more 
readily  available  than  geoacoustic  measurements  and  BL  measure¬ 
ments.  Thus,  the  first  prerequisite  for  Biot/Stoll  modeling  to 
improve  capability  in  BL  prediction  has  been  met. 

Examples  can  be  cited  for  significantly  reduced  variances  in 
some  provinces  compared  to  ZZ.  Especially  dramatic  are  the  cases 
of  provinces  AA,  AB  and  AC.  All  show  significantly  greater  pre¬ 
ponderance  of  low  BL  conditions  and  significantly  lower  variance 
of  BL  values  than  does  ZZ.  Not  only  are  variances  reduced  in 
some  provinces,  but  means  are  significantly  different  in  others 
as  shown  in  Figures  5-2,  5-4,  and  5-9.  Thus,  the  second  pre¬ 
requisite  for  Biot/Stoll  modeling  to  improve  capability  in  BL 
prediction  has  also  been  met. 

An  insufficient  amount  of  compressional  speed  data  has  pre¬ 
cluded  our  formally  estimating  the  probability  that  a  sediment 
was  incorrectly  assigned  to  a  province.  In  light  of  our  other 
findings,  however,  we  do  not  need  such  an  estimate  to  show  sub¬ 
stantial  improvement  due  to  provincing. 

5.4  Conclusions  and  Recommendations 

We  conclude  that  PHYSED  modeling,  when  combined  with 
provincing  based  on  water  depth  and  qualitative  size  class, 
represents  a  real  improvement  in  our  capability  to  predict  Bottom 
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Loss  (BL)  on  the  continental  terrace.  This  conclusion  is  based 
both  on  reduced  variances  in  some  provinces  and  on  relocated 
means  in  some  provinces. 

We  also  conclude  that  the  provincing  done  here  provides  a 
better  separation  from  combined  sediments  (province  ZZ)  for  AA, 
AB  and  AC  (shelf  sands,  shelf  dirty  sands,  and  shelf  sandy  muds) 
than  it  does  for  the  other  provinces.  This  is  based  upon  the 
significant  differences  found  in  the  proportion  of  low  BL,  low  BL 
means,  and  low  BL  variances. 

We  recommend  that  the  sediments  in  the  other  provinces 
should  be  divided  in  different  ways  than  done  here  in  an  attempt 
to  discover  a  provincing  scheme  that  reduces  the  variance  under 
high  BL  conditions.  Separation  by  ocean  area  is  suggested  by  our 
void  ratio  analysis  in  Section  2.2.  Separation  by  void  ratio 
ranges  directly  is  also  reasonable  but  this  procedure  depends  on 
engineering  data  which  is  less  readily  available  than  ocean  area. 

Before  pursuing  further  provincing  schemes,  we  recommend 
increasing  the  amount  of  data  in  PHYPROSE.  Province  BA  seems  to 
have  been  affected  by  unrepresentative  high  grain  densities. 
Also,  the  near  agreement  of  the  BL  distribution  belonging  to 
province  ZZ  with  province  BE,  under  high  BL  conditions,  points  to 
a  high  percentage  of  silts  in  the  PHYPROSE  data  base. 
Continental  terrace  sediments  as  a  whole  were  expected  to  be 
closer  to  sands.  The  high  silt  content  in  ZZ  may  be  an  artifact 
of  using  data  from  all  depths  of  a  core,  which  results  in  more 
observations  per  silty  site  than  observations  per  sandy  site. 
Another  reason  for  this  possible  error  is  that  the  data  selected 
may  not  be  a  random  sample  of  real  terrace  sediments  but  are,  for 
some  reason,  biased. 

Given  these  provincing  results  and  the  current  interest  in 
developing  geoacoustic  or  bottom  loss  models  for  operational  use 
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in  shallow  water  areas,  we  strongly  recommend  the  continued 
development  of  the  physical  sediment  approach.  In  Section  1.1  we 
discussed  the  stages  of  development  of  this  hybrid  model.  We 
have  shown  significant  progress  through  six  stages  of  the  seven 
stage  development  plan  displayed  as  Table  1-4.  Progress  on  Stage 
VI,  Predictive  Model  Evaluation  had  been  shown  in  this  effort. 
However,  the  scarcity  of  acoustic  measurements  precluded  a 
credible  model-to-data  evaluation.  Effort  should  be  directed 
toward  such  an  evaluation.  Once  we  are  confident  of  the  model, 
then  validation  efforts  (Stage  VII)  will  be  warranted.  We  exoect 
to  continue  such  efforts  under  Navy  research  and  development. 
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FORTRAN  CODE  for  major  programs  used  in 
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PROGRAM  SAUCE 


Program,  based  on  program  Tomato,  to  select  data  from 
provinces  defined  by: 

1. )  Water  Depth 

2. )  Sediment  Size  Type 

3 .  )  Ocean  Area 


Provinces  will  be  indicated  using  four  characters  as  follows: 
1st  character  indicates  Water  Depth. 

A  =>  0  to  200  meters 
B  *  201  to  2000  meters 
C  =  2001  meters  and  greater 
Z  =  All  depths  combined 


2nd  character  indicates  sediment  size  type. 

A  =  Sand  (Qual.  Size  code  1) 

Sand  with  finer  fractions  (codes  2,  3,  and  4) 
Sandy  finer  fractions  (codes  5,  6,  and  7) 

Silt  (code  8) 

Clays  (codes  9  and  10) 

All  sizes  combined 


B 

C 

D 

E 

Z 


3rd  &  4th  characters  indicate  Ocean  Area  codes  (see  Phyprose 
Report) . 

ZZ  =  All  Ocean  Areas  Combined 


Output  files  willbe  named  with  Province  Name  and  Parameter 
codes  as  follows: 

****VOID  *  Void  ratios  in  province 
****SGGR  =  Grain  densities  in  province 
****SZRO  =  So,  mean  specific  surfaces  in 
province 


Where  ****  indicates  the  province  name 


This  naming  convention  can  easily  be  expanded  to  include 
other  provincing  criteria,  as  needed. 


LOGICAL* 2  WORK 

character* 14  FNAM1,  FNAM2,  FNAM3 
character *13  FNAM4 

character* 1  IALPH(6,2) , I CHAR (2) 
character* 10  STRING, NCRUIS, NS AMP 
common/Al/ARRl , OUT1 
common/ A2/ARR2 , 0UT2 
common/A3/ARR3 , 0UT3 

dimension  ARR1(5000,2) ,OUT1(100,2) , ARR2 (5000, 2) ,OUT2(100,2) 
+  ARR3 (9500 , 2 ) ,OUT3(100,2) ,GSPD(16) , TEMPV ( 100 , 2 ) , 

+  TEMPZ (100, 2) , NACV (100) ,NACZ (100) 


data  IALPH/'A'  /B'/C'/D'/E'/Z'/a'/b'/c'/d'/e'/z'/ 


7001  format  ( 
+ 

+ 

+ 
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ENTER  WATER  DEPTH:’/ 

A  =  0  to  200  meters'/ 

B  =  201  to  2000  meters'/ 

C  =  2001  to  and  greater'/ 
Z  =  All  depths') 
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7002 


7003 


6000 

6001 

6002 

6010 

6011 

6012 

6013 

6030 

6031 


format 


format 


format 

format 

format 

format 

format 

format 

format 

format 

format 

format 

format 

format 

format 

format 

format 

format 


ENTER  SEDIMENT  SIZE  TYPE:'/ 

A  =  Sand  (Code  1) •/ 

B  =  Sand  &  Finers  (Codes  2,3,&  4) ', 
C  =  Sandy  Finers  (Codes  5,6,&  7)'/ 
D  =  Silt  (Code  8) '/ 

E  =  Clays  (Codes  9  &  10) '/ 

Z  =  All  Sizes  (Codes  1  thru  10)') 
ENTER  OCEAN  AREA: '/ 

Two  digit  Ocean  Area  code  or'/ 
enter  ZZ  for  all  ocean  areas'/ 
(Codes  00  thru  99)') 

WATER  DEPTH : ' ) 

SED.  SIZE  TYPE:') 

OCEAN  AREA: • ) 

3x, '  0  to  200  meters  ') 

3x, '  201  to  2000  meters  ') 

3x, '  2001  meters  &  greater  ') 

3x, '  All  Depths  ') 

5x,  1  Sand  (Code  1)  ') 

:  I  f.  T7  -i  n  ->  r  a  \  • 


(13x, ' 

( 13x, ' 

( 13x, ' 

(13X,  ' 

(15x, '  sand  (code  l)  ') 

( 15x, 1  Sand  &  Finers  (Codes  2,3,&4)  ') 
( 15x, '  Sandy  Finers  (Codes  5,6,&7)  ') 

( 15x, '  Silt  (Code  8)  ') 

( 15x, '  Clays  (Codes  9&10)  ') 

(15x, '  All  Sizes  (Codes  1  thru  10)  ') 

(13x,'  All  areas  (Codes  00  thru  99)  ') 

(13x,i3) 

/  /  /  \ 


RK  =  .FALSE. 

write  (*,6099) 
write  (*,7001) 
read  (*,6101)  ICHAR(l) 
format  (al) 

do  50  I  =  1,  6 
do  50  J  =  1,  2 

if  (ICHAR(l)  .eq.  IALPH(I , J) )  goto  60 
continue 
write  (*,6102) 

format  (/'  *?*  NOT  ACCEPTABLE  RESPONSE.  TRY  AGAIN.'/) 
goto  30 

if  (I  .eq.  4  .or.  I  .eq.  5)  goto  55 
LOOPD  =  I 

write  (*,6099) 
write  (*,7002) 
read  (*,6101)  ICHAR(l) 

do  100  I  =  1,  6 
do  100  J  =  1,  2 

if  (ICHAR(l)  .eq.  IALPH(I,J))  goto  110 
continue 
write  (*,6102) 
goto  80 
LOOPS  =  I 

write  (*,6099) 
write  (*,7003) 

read  (*,6104)  ICHAR(l) , ICHAR(2) 


SAUCE 


6104  format  (2al) 

if  ((ICHAR(l)  .eq.  IALPH(6,1)  .or.  ICHAR(l)  .eq.  IALPH(6,2) 


rS' 

+ 

.  and. 

& 

+ 

(ICHAR (2)  .eq.  IALPH(6,1) 

+ 

then 

LOOPA  =  0 

else 

if  (ICHAR(2)  .eq.  •  ')  then 
ICHAR(2)  =  I CHAR ( 1) 
if  ( I CHAR ( 2 )  .eq.  '  •)  then 
I CHAR ( 1)  =  'Z' 

ICHAR(2)  =  'Z' 
else 

I CHAR ( 1)  =  '  1 
end  if 
end  if 

write  (string, 6104)  ICHAR(i) , ICHAR(2) 
read  (STRING, 6105 , ERR=135)  10 ARE A 

6105  format  (i2) 

LOO PA  =  1 

end  if 
goto  140 

135  write  (*,6102) 
goto  130 

140  continue 

Generate  output  file  names  in  FNAM1,  FNAM2,  FNAM3 

if  (LOOPA  .eq.  0)  then 
ICHAR(l)  =  *Z' 

I CHAR ( 2 )  =  'Z' 
end  if 

write  (FNAM1, 6106) IALPH(LOOPD, 1) , IALPH ( LOOPS , 1 ) , 

+  I CHAR  ( 1 ) , I CHAR ( 2 ) 

6106  format  (' E: ', 4al, 'VOID. DAT ' ) 

write  (FNAM2, 6107) IALPH (LOOPD,l) , IALPH (LOOPS , 1) , 

+  I CHAR ( 1 ) , ICHAR ( 2 ) 

6107  format  ( 'E: ' ,4al, 'SGGR.DAT' ) 

write  ( FNAM3 , 6108 ) IALPH ( LOOPD, 1) , IALPH (LOOPS , 1) , 

+  ICHAR(l) ,ICHAR(2) 

6108  format  ( 'E: ' ,4al, 'SZR0.DAT' ) 

WRITE  ( FNAM4 ,6111)  IALPH (LOOPD, 1) , IALPH (LOOPS , 1) , 

+  ICHAR(l) ,ICHAR(2) 

6111  FORMAT  ( ' E : ' , 4 A1 , ' COV . DAT ' ) 

List  results  for  operator. 

write  (*,6109) 

6109  format  (///'  YOU  CHOSE  A  PROVINCE  WITH  ') 
write  (*,6000) 

if  (LOOPD  .eq.  1)  then 
write  (*,6010) 
else  if  (LOOPD  .eq.  2)  then 
write  (*,6011) 
else  if  (LOOPD  .eq.  3)  then 
write  (*,6012) 
else 
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write  (*,6013) 
end  if 

write  (*,6001) 
if  (LOOPS  .eq.  1)  then 
write  (*,6030) 
else  if  (LOOPS  .eq.  2)  then 
write  (*,6031) 
else  if  (LOOPS  .eq.  3)  then 
write  (*,6032) 
else  if  (LOOPS  .eq.  4)  then 
write  (*,6033) 
else  if  (LOOPS  .eq.  5)  then 
write  (*,6034) 
else 

write  (*,6035) 
end  if 

write  (*,6002) 
if  (LOOPA  .eq.  0)  then 
write  (*,6050) 
else 

write  (*,6051)  IOAREA 
end  if 

write  (*,6110)  FNAM1,  FNAM2 ,  FNAM3 
format  (/ '  FOR  THIS  PROVINCE'/ 

'  VOID  RATIOS  ARE  IN  FILE  ' ,A12/ 

'  GRAIN  DENSITIES  ARE  IN  FILE  ' ,A12/ 

'  &  SPECIFIC  SURFACE  MEANS  ARE  IN  FILE  ’,A12) 

Open  output  files 

open  (60,  file=FNAMl,  status* ' new 1 , 

access* ' sequential ' ,  form* ' formatted ' ) 
open  (70,  file=FNAM2,  status* ' new ' , 

access* ' sequential ' ,  form* 1  formatted ' ) 
open  (80,  file=FNAM3,  status* ' new' , 

access* ' sequential ' ,  form*' formatted ' ) 

OPEN  (90 , FILE=FNAM4 ,  STATUS* ' NEW ' , 

-  ACCESS  =  'SEQUENTIAL',  FORM  =  'FORMATTED') 

Open  input  files 

open  (10,  file  =  'HEADER. DAT' ,  status  =  'OLD', 
i-  access  =  'SEQUENTIAL',  form  =  'UNFORMATTED') 

open  (20,  file  =  'C0MM0NP.DAT',  status  =  'OLD', 
t-  access*  '  DIRECT '  ,  form*  '  UNFORMATTED '  ,  recl=44 ) 

open  (30,  file  =  'SIZEDIS.DAT',  status  =  'OLD', 
f  access  *  'DIRECT',  form  =  'UNFORMATTED',  reel  =  80) 

open  (40,  file  =  'SIZECLS.DAT',  status  =  'OLD', 
f  access  =  'DIRECT',  form  =  'UNFORMATTED',  reel  =  36) 

open  (50,  file  =  ' C : QUALS I Z . DAT ' ,  status  =  'OLD', 

+  access  =  'DIRECT',  form  =  'UNFORMATTED',  reel  =  20) 

NVOID  =  0 
NSGGR  =  0 
NSZRO  =  0 
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KOUNTV  =  0 
KOUNTZ  =  0 

200  read  (10,  end  =  500)  NAC, 

+  INSTS , IPLAT , NCRUIS , NSAMP , LAT , DLAT , LON , DLON , 

+  IOCEAN , IDAY , MONTH, I YEAR, MEAST , COREL, WDEPTH , 

+  IWDM, IPHYSG, ISURST, ICOMP, IENGP, IMINP, IQGS , 

+  IGSC, IGSD, IACOP, IFLUD 

201  format  (i6,  i4 ,  lx,  i.4 ,  lx,  alO ,  alO ,  i3 ,  f  5 . 4  ,  i4 ,  f 5 . 4  , 

+  Ix,4i2,i3,f7.2,f6.1,2x,3i2,5x,8i6) 

C  Reject  if  out  of  area,  unless  all  areas  are  accepted. 

if  (LOOPA  .eq.  0)  goto  220 
if  (IOCEAN  .ne.  IO ARE A)  goto  200 

C  Reject  if  out  of  water  depth,  unless  all  depths  are  accepted. 

220  if  (LOOPD  .eq.  6)  goto  230 

if  (WDEPTH  .It.  0.0)  goto  200 
if  (LOOPD  .eq.  1)  then 

if  (WDEPTH  .gt.  200.0)  goto  200 
else  if  (LOOPD  .eq.  2)  then 

if  (WDEPTH  .le.  200.0  .or.  WDEPTH  .gt.  2000.0)  goto  200 
else  if  (LOOPD  .eq.  3)  then 

if  (WDEPTH  .le.  2000.0)  goto  200 
end  if 

C  Reject  if  no  Sed.  Size  info,  unless  all  types  are  accepted. 

C  NOTE:  Must  do  final  reject  on  size  type  at  each  depth  in  core. 

230  if  (LOOPS  .eq.  6)  goto  240 

if  (IQGS  .eq.  0  .and.  IGSC  .eq.  0  .and.  IGSD  .eq.  0)  goto  200 

C  Load  VOID  RATIO  and  SGGR  values  into  ARRl  and  ARR2 . 

240  if  (ICOMP  .eq.  0)  goto  400 

300  read  (20,  rec  =  ICOMP,  err=500)  NACCES,DTOP,DBOT, 

+  WETWT , SGGR , WATCON , VO I D , S VO I D , POROS , DMG , PERM 


301  format  (i6, 2f7 . 2 , 2f 4 . 2 , f 7 . 1 , f 6 . 3 , f 7 . 3 , f 4 . 3 , f7 . 4 , eio .  3) 

if  (NAC  .ne.  NACCES)  goto  400 
if  (LOOPS  .eq.  6)  goto  340 
IGSD1  =  IGSD 
IGSC1  =  IGSC 
IQGS1  =  IQGS 

call  TYPES (NAC , IGSD1 , IGSC1 , IQGS 1 , ICLASS , DTOP , DBOT) 
if  (ICLASS  .eq.  12)  goto  310 
if  (LOOPS  .eq.  1)  then 

if  (ICLASS  .ne.  1)  goto  310 
else  if  (LOOPS  .eq.  2)  then 

if  (ICLASS  .ne.  2  .and.  ICLASS  .ne.  3  .and. 

+  ICLASS  .ne.  4)  goto  310 

else  if  (LOOPS  .eq.  3)  then 

if  (ICLASS  .ne.  5  .and.  ICLASS  .ne.  6  .and. 

+  ICLASS  .ne.  7)  goto  310 

else  if  (LOOPS  .eq.  4)  then 
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if  (ICLASS  . ne.  8)  goto  310 
else  if  (LOOPS  .eq.  5)  then 

if  (ICLASS  .ne.  9  .and.  ICLASS  .ne.  10)  goto  310 
end  if 
goto  330 

310  ICOMP  =  ICOMP  +  1 
goto  300 
330  continue 

C  Present  class  is  correct,  load  VOID  RATIO  &  SGGR  if  present. 

C  1st,  VOID  RATIO 

340  if  (VOID  .ne.  -9.0)  then 
NVOID  =  NVOID  +  1 
ARR1 (NVOID, 1)  =0.0 
ARR1 (NVOID, 2 )  =  VOID 
KOUNTV  =  KOUNTV  +  1 
TEMPV ( KOUNTV , 1 )  =  DTOP 
TEMPV (KOUNTV, 2)  =  VOID 
NACV (KOUNTV)  =  NACCES 
else  if  (SVOID  .ne.  -9.0)  then 
NVOID  =  NVOID  +  1 
ARR1 (NVOID, 1)  =0.0 
ARR1 (NVOID, 2 )  =  SVOID 
KOUNTV  =  KOUNTV  +  1 
TEMPV ( KOUNTV, 1)  =  DTOP 
TEMPV (KOUNTV, 2)  =  SVOID 
NACV (KOUNTV)  =  NACCES 
else  if  (POROS  .gt.  0.0)  then 
NVOID  =  NVOID  +  1 
ARR1 (NVOID, 1)  =0.0 

ARR1 (NVOID, 2 )  =  POROS  /  (1.0  -  POROS) 

KOUNTV  =  KOUNTV  +  1 
TEMPV ( KOUNTV , 1 )  =  DTOP 

TEMPV (KOUNTV, 2)  =  POROS  /  (1.0  -  POROS) 

NACV (KOUNTV)  =  NACCES 
end  if 

Now,  SGGR 

350  if  (SGGR  .gt.  0.0)  then 

write  (*,*)  'Nacces  =  ',nacces 

write  (*,*)  ' Dtop  =  ' ,dtop 

write  (*,*)  ' Dbot  =  ',dbot 

NSGGR  =  NSGGR  +  1 

ARR2 (NSGGR, 1)  =0.0 
ARR2 (NSGGR, 2 )  =  SGGR 
end  if 

ICOMP  =  ICOMP  +  1 
goto  300 

Now  read  FULL  GRAIN  SIZE  DISTRIBUTION  records  to  compute  SZRO 

400  if  (IGSD  .eq.  0)  then 
KOUNTV  =  0 
GOTO  200 

end  if 

401  read  (30,  rec  =  IGSD,  err  =  200)  NA, 

+  DTOP, DBOT, GSPD(l) ,GSPD(2) ,GSPD(3) ,GSPD(4) , 


SAUCE 


GSPD (5) , GSPD (6) ,GSPD(7) fGSPD(8) ,GSPD(9) ,GSPD(10) , 
GSPD(ll) , GSPD (12) , GSPD ( 13 ) ,GSPD(14) ,GSPD(15) , 
GSPD (16) ,DM 


402  format  (16, 2F7 . 2 , 16F7 . 3 , F7 . 4) 
if  (NA  .ne.  NAC)  THEN 

IF  (KOUNTV  .GT.  0  .AND.  KOUNTZ  .GT.  0)  THEN 
1  =  1 
J  =  1 

405  IF  (TEMPV(I , 1)  . EQ.  TEMPZ(J,1))  THEN 

WRITE  (90,6112)  NACV(I) , TEMPV(I , 2 ) , NACZ (I) , TEMPZ ( J, 2) 
6112  FORMAT  (15 , IX, F7 . 3 , IX, 15 , G13 . 5) 

IF  (I  .EQ.  KOUNTV  .AND.  J  .EQ.  KOUNTZ)  GOTO  408 
IF  (I  .NE.  KOUNTV)  1=1+1 
IF  (J  .NE.  KOUNTZ)  J  =  J  +  1 
ELSE  IF  (TEMPV (1,1)  .LT.  TEMPZ (J,l))  THEN 
IF  (I  .EQ.  KOUNTV)  GOTO  408 
1  =  1  +  1 

ELSE  IF  (TEMPV (1,1)  .GT.  TEMPZ (J,l))  THEN 
IF  (J  .EQ.  KOUNTZ)  GOTO  408 
J  =  J  +  1 
END  IF 
GOTO  405 
408  CONTINUE 
END  IF 
KOUNTV  =  0 
KOUNTZ  =  0 
GOTO  200 
END  IF 

if  (LOOPS  .eq.  6)  goto  440 
IGSD1  =  IGSD 
IGSC1  =  IGSC 
IQGS1  =  IQGS 

call  TYPES (NAC, IGSD1 , IGSC1 , IQGS1 , ICLASS , DTOP, DBOT) 
if  (ICLASS  .eq.  12)  goto  410 
if  (LOOPS  .eq.  1)  then 

if  (ICLASS  .ne.  1)  goto  410 
else  if  (LOOPS  .eq.  2)  then 

if  (ICLASS  .ne.  2  .and.  ICLASS  .ne.  3  .and. 


if  (ICLASS  .ne.  1)  goto  410 
else  if  (LOOPS  .eq.  2)  then 

if  (ICLASS  .ne.  2  .and.  ICU 
ICLASS  .ne.  4)  goto  410 
else  if  (LOOPS  .eq.  3)  then 


,  ne .  3  . and . 


if  (ICLASS  .ne.  5  .and.  ICLASS  .ne.  6  .and. 

ICLASS  .ne.  7)  goto  410 
e  if  (LOOPS  .eq.  4)  then 
if  (ICLASS  .ne.  8)  goto  410 
e  if  (LOOPS  .eq.  5)  then 

if  (ICLASS  .ne.  9  .and.  ICLASS  .ne.  10)  goto  410 


else  if  (LOOPS 
if  (ICLASS  . 
else  if  (LOOPS 
if  (ICLASS  . 
end  if 
goto  430 
IGSD  =  IGSD  +  1 
goto  401 
continue 


Present  class  is  correct,  compute  SZRO 

SZRO  =0.0 
do  450  1=1,  16 

PHI  =  (float (I)  +  0.5)  -  6.0 
if  (nac  .eq.  34)  write  (*,*)  i,phi 


DIAM  =  (2.0  **  (-  PHI))  /  10.0 
if  (nac  .eq.  34)  write  (*,*)  diam 
SO  *  6.0  /  DIAM 

if  (nac  .eq.  34)  write  (*,*)  s0,gspd(i) 
if  (GSPD(I)  .gt.  0.0)  then 

SZRO  =  SZRO  +  SO  *  GSPD(I)  /  100.0 
if  (nac  .eq.  34)  write  (*,*)  szro 
end  if 

450  continue 

if  (nac  .eq.  34)  write  (*,*)  szro 
NSZRO  =  NSZRO  +  1 
ARR3 (NSZRO, 1)  =0.0 
ARR3 (NSZRO, 2 )  =  SZRO 

KOUNTZ  =  KOUNTZ  +  1 

TEMPZ (KOUNTZ , 1)  =  DTOP 

TEMPZ (KOUNTZ, 2)  =  SZRO 

NACZ (KOUNTZ)  =  NA 
IGSD  =  IGSD  +  1 
goto  401 

C  Arrays  are  full  of  data,  now  compute  cumulative 

C  distribution  functions  and  write  files. 

500  CONTINUE 

CLOSE  (90) 

C  Void  Ratio 

WRITE  (*,*)  'NVOID  =  ' ,NVOID 

call  CDF (ARR1, 5000, NVOID, OUTl,Nl) 
write  (60,601)  Nl, (0UT1(J,1) ,OUTl(J,2) ,J=1,N1) 
601  format  (lx,  i3 , 100  ( f7 . 4 ,  f7 . 3) ) 

CLOSE  (60) 

C  Grain  Density 

WRITE  (*,*)  ' NSGGR  =  1 , NSGGR 

call  CDF (ARR2, 5000, NSGGR, OUT2,N2) 

write  (70,601)  N2 , (0UT2 (J, 1) ,OUT2 (J, 2) , J=1,N2) 

CLOSE  (70) 

C  Mean  Specific  Surface 

WRITE  (*,*)  'NSZRO  =  ',  NSZRO 

call  CDF (ARR3, 9500, NSZRO, OUT3,N3) 
write  (80,801)  N3 , (0UT3 ( J , 1) , OUT3 ( J, 2 ) , J=1 , N3 ) 
801  format  (lx, i3 , 100 (f7. 4 , ell. 4) ) 

CLOSE  (80) 

Close  all  files 

close  (10) 
close  (20) 
close  (30) 
close  (40) 
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close  (50) 

stop 

K- 

end 

C* ************************************************************ 

SUBROUTINE  CDF(ARR, ISIZE,N, OUT, NOUT) 

C  Takes  data  in  the  array  ARR,  computes  a  cumulative 

C  distribution  function  (CDF)  and  returns  the  CDF  in 

C  the  array  OUT.  ARR  has  N  values  (up  to  9500) ,  but 

C  OUT  will  contain  NOUT  values  (no  more  than  100) . 

dimension  ARR (ISIZE, 2) ,  OUT(100,2) 

C  "Bubble”  sort  routine  to  rank  order  the  measurements 

C  from  1=1  (being  the  largest)  to  I  =  N  (being  the  smallest) . 

if  (N  .eq.  1)  then 
NOUT  =  1 
OUT (1,1)  =  0.0 
OUT (1,2)  =  ARR(1, 2) 

GOTO  150 
end  if 
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if  (ARR (2,2)  .gt.  ARR(1,2))  then 
A  =  ARR (2,2) 

ARR (2,2)  =  ARR(1,2) 

ARR (1,2)  =  A 
end  if 


if  (N  .eq.  2)  then 
NOUT  =  2 
OUT (1,1)  =0.0 
OUT (2,1)  =  1.0 
OUT (1,2)  =  ARR(2 , 2) 
OUT (2,2)  =  ARR( 1 , 2 ) 
GOTO  150 
end  if 


C  WRITE  (*,*)  '  ENTER  110  LOOP* 

do  110  I  =  3,  N 

if  (ARR(I, 2)  .gt.  ARR(I-l, 2) )  then 
A  =  ARR ( 1-1 , 2 ) 

ARR(I-1, 2)  =  ARR (1,2) 

ARR (1,2)  =  A 
K  =  1 

70  if  (ARR(I-K, 2)  .gt.  ARR ( I- (K+l) , 2 ) )  then 

A  =  ARR (I- (K+l) ,2) 

ARR(I-(K+1) ,2)  =  ARR ( I-K, 2 ) 

ARR ( I-K, 2 )  =  A 
K  =  K  +  1 

if  (I-K  .eq.  1)  goto  110 
goto  70 
end  if 
end  if 

110  continue 

C  WRITE  (*,*)  '  EXIT  110  LOOP' 
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Now  compute  the  probabilities  for  each  observed  value. 

The  array  ARR(I, J=l, 2)  is  now  ordered  from  1=1  (Largest) 
to  I=N  (Smallest) .  Duplicate  values  presumbably  are 
located  adjacent  to  one  another. 

WRITE  (*,*)  '  ENTER  120  LOOP* 

ARR(1,1)  =  1.0000 
do  120  I  -  2,  N 

ARR(I, 1)  =  FLOAT ( N  -  (I  -  1))  /  FLOAT ( N ) 

120  continue 

WRITE  (*,*)  '  EXIT  120  LOOP* 

Set  to  zero  the  probability  value  for  all  but  the  first 
observation  in  a  string  of  duplicates. 

WRITE  (*,*)  '  ENTER  130  LOOP1 

KICK  =  0 

do  130  I  =  1,  N  -  1 

if  (ARR(I , 2)  .eq.  ARR(I+1,2))  then 
ARR(I+1, 1)  =  0.0 
KICK  =  KICK  +  1 
end  if 

130  continue 

WRITE  (*,*)  '  EXIT  130  LOOP* 

If  more  than  99  points  still  exist,  then  set  the  probabilities 
of  the  values,  which  are  cloest  together,  to  zero. 

(Don't  drop  the  first  or  last  points.) 

WRITE  (*,*)  '  ENTER  132  to  139  LOOP' 

132  DMIN  =  ARR (1,2)  -  ARR(N, 2) 

WRITE  (*,*)  ' DMIN= ' , DMIN , '  ARR( 1 , 2 ) = • , ARR ( 1 , 2 ) , 

+  '  ARR (N , 2 ) = ' , ARR (N , 2 ) 

WRITE  (*,*)  '  N= ' , N , '  KICK= ' , KICK 

if  ((N  -  KICK)  .le.  99)  goto  139 
1  =  1 

133  J  =  1 

134  CONTINUE 

WRITE  (*,*)'  I-', I,'  J=',J 

if  (I+J  .ge.  N)  goto  136 
WRITE  (*,*)  '  ARR ( I + J , 1 ) = ' , ARR ( I + J , 1 ) 

if  (ARR(I+J,1)  .ne.  0)  goto  135 
J  =  J  +  1 

WRITE  (*,*)  'J=',J 

goto  134 

135  DIF  =  ARR(I , 2)  -  ARR(I+J,2) 

V— TE  (*,*)  '  DIF=  '  ,  DIF,  '  ARR (1,2)  =  '  ,  ARR (1,2)  , 

'  ARR (I+J, 2) =' , ARR (I+J, 2) 

WRITE  (*,*)  '  DIF*', DIF,'  DMIN= ' , DMIN 

if  (DIF  .It.  DMIN)  then 
IMIN  =  I  +  J 
DMIN  =  DIF 
end  if 
I  =  I  +  J 
goto  133 

136  ARR(IMIN, 1)  =  0.0 

KICK  =  KICK  +  1 
write  (*,*)  kick 
goto  132 
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139  continue 

C  WRITE  (*,*)  '  EXIT  132  to  139  LOOP* 

Next  condense  the  array  to  "plot"  only  those  points  with 
non-zero  frequency.  Condense  from  smallest  value  to 
largest  value. 

WRITE  (*,*)  '  ENTER  140  LOOP* 

J  -  1 

do  140  K  =  1,  N 
I  =  N  +  1  -  K 

if  (ARR(I, 1)  .ne.  0.0)  then 
J  *  J  +  1 

OUT(J,l)  =  ARR(I , 1) 

OUT ( J ,  2 )  =  ARR(I, 2) 
end  if 

140  continue 
NOUT  =  J 

WRITE  (*,*)  '  EXIT  140  LOOP' 

Fill  OUT (1,2)  with  a  value  extrapolated  from  overall 
distribution  slope  to  "zero"  probability. 

SLOPE  =  (OUT (NOUT, 2)  -  OUT(2,2))  /  (1.0  -  OUT(2,l)) 

OUT (1,2)  =  OUT ( NOUT , 2 )  -  SLOPE 
OUT (1,1)  =  0.0 

150  return 
end 

C* ************************************************************ 


subroutinetypes(nac, igsd, igsc, iqgs, iclass,dtop,dbot) 
dimension  gspd(16) 

10  format  ('  ',i6,'  ',i6,'  ',i6,'  ',i6,'  ',i2,'  ',f7.2,'  ' , f 7 .2) 

1000  PCOURSE  =0.0 

PSAND  =0.0 
PSILT  =0.0 
PCLAY  =0.0 
PFINE  =0.0 

if  (IGSD  .gt.  0)  then 

1400  read  (30,  rec  =  IGSD,  err  =  1406)  NA, 

+  DT0P1 , DBOT1 , GSPD ( 1)  ,GSPD(2)  ,GSPD(3)  ,GSPD(4)  , 

+  GSPD (5) ,GSPD(6) ,GSPD(7) ,GSPD(8) ,GSPD(9) ,GSPD(10) 

+  GSPD (11) , GSPD ( 12 ) , GSPD ( 13 ) ,GSPD(14) ,GSPD(15) , 

+  GSPD (16), DM 

1401  format  (16 , 2F7 . 2 , 16F7 . 3 , F7 . 4 ) 

if  (NA  .ne.  NAC)  goto  1406 

if  (DTOP  .ne.  DTOP1)  then 
if  (DBOT1  .le.  DTOP  .or.  DTOP1  .ge.  DBOT)  then 
IGSD  =  IGSD+1 


GOTO  1400 
end  if 
end  if 

do  1402  K  =  1,  4 

if  (GSPD(K)  .ge.  0)  PCOURSE  =  PCOURSE  +  GSPD(K) 
continue 

do  1403  K  =  5,  9 

if  (GSPD(K)  .ge.  0)  PSAND  =  PSAND  +  GSPD(K) 
continue 

do  1404  K  =  10,  13 

if  (GSPD(K)  .ge.  0)  PSILT  =  PSILT  +  GSPD(K) 
continue 

do  1405  K  =  14,  16 

if  (GSPD(K)  .ge.  0)  PCLAY  =  PCLAY  +  GSPD(K) 
continue 

goto  1800 

end  if 

if  (IGSC  .gt.  0)  then 

read  (40,  rec  =  IGSC,  err  =  1502)  NA, 

DTOP2 , DBOT2 , PCOURSE , PSAND , PSILT , PCLAY , PMUD , DD 

format  (16,  2F7.2,  5F7.3,  F7.4) 

if  (NA  .ne.  NAC)  goto  1502 

if  (DTOP  .ne.  DT02)  then 
if  ( DBOT2  .le.  DTOP  .or.  DTOP2  .ge.  DBOT)  then 
IGSC  =  IGSC  +  1 
goto  1500 
end  if 
end  if 

goto  1800 

end  if 

if  (IQGS  .gt.  0)  then 

read  (50,  rec  =  IQGS,  err  =  1700)  NA, 

DTOP3 ,  DB0T3 ,  DMG,  IQSC 

format  (i6,  f 7 . 2 ,  f7.2,  f7.4,  i2) 

if  (NA  .ne.  NAC)  goto  1700 

if  (DTOP  .ne.  DTOP3 )  then 

if  (DB0T3  .le.  DTOP  .or.  DTOP3  .ge.  DBOT)  then 
IQGS  =  IQGS  +  1 
goto  1600 
end  if 
end  if 
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ICLASS  =  IQSC 
goto  1900 

end  if 

ICLASS  -  12 
goto  1900 

PFINE  =  PS AND  +  PSILT  +  PC LAY 

if  (PCOURSE  .gt.  50.0)  then 
ICLASS  =  11 
goto  1900 
end  if 

if  (PFINE  .It.  0.001)  then 
ICLASS  =  12 
goto  1900 
end  if 

PI  =  PSAND  /  PFINE 

if  (PI  .ge.  0.75)  then 
ICLASS  =  1 
goto  1900 
end  if 

if  (PI  .It.  .75  .and.  PI  .ge.  .5)  then 
ICLASS  =  3 

P2  -  PSILT  /  (PSILT  +  PCLAY) 

if  (P2  .It.  0.33)  ICLASS  =  4 
if  (P2  .gt.  0.66)  ICLASS  =  2 

goto  1900 
end  if 

if  (PI  .It.  0.5  .and.  PI  .ge.  0.25)  then 
ICLASS  =  6 

P2  =  PSILT  /  (PSILT  +  PCLAY) 

if  (P2  .It.  0.33)  ICLASS  =  7 

if  (P2  .gt.  0.66)  ICLASS  =  5 

goto  1900 
end  if 

if  (PI  .It.  0.25)  then 
ICLASS  =  9 

P2  *  PSILT  /  (PSILT  +  PCLAY) 

if  (P2  .It.  0.33)  ICLASS  =  10 

if  (P2  .gt.  0.66)  ICLASS  =  8 

end  if 

continue 

return 

end 
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PROGRAM  MCPHYS 


A  NEW  VERSION  OF  PHYSED  ALTERED  FOR  THE 
ALLOWANCE  OF  MONTE  CARLO  TESTING  OF  THE  VARIABILITY 
OF  THE  BIOT-STOLL  THEORY  IN  PHYSED.  STATISTICS 
FOR  THE  TESTING  ARE  DERIVED  FROM  THE  PHYPROSE 
DATA  BASE  AND  APPLY  TO  THE  FOLLOWING  VARIABLES: 

RHOR  (  GRAIN  DENSITY  ) 

VOID  (  VOID  RATIO  ) 

SO  (  MEAN  SPECIFIC  SURFACE  ) 

DECBS  (  SHEAR  LOG  DECREMENT  ) 

DECBE  (  COMPRESS IONAL  LOG  DECREMENT  ) 

PRATIO  (  POISSON'S  RATIO  ) 


COMMON/ CONSTS/BLKRR , RHOF , BLKFR , ETA , ALPHA , HOVEMK , 

+  GC1,GC2,GC3,VDH 

COMMON/VDEPTH/EVERY (20) 

COMMON/ VFREQ/ FREQ ( 16 ) 

COMMON/ DISTVALS/RHOR, VOID, SO , DECBS , DECBE, PRATIO 
COMMON / FCTVLS / BETA , PERM, PSIZE 
COMMON/DEPLOOP/EPRESS , G , GPRI , E , EPRI 

COMMON/BCOEF/CBLKR, CBLKF , CBLKB , CSHRB , CDEN , CFAC1 , CH , CC , CM 

COMMON/FVARS/OMEGA, EKAPPA, EM 

COMMON/ COMKLV/CFAC , CFAC2 , CFAC3 , CFAC4 , CF 

COMMON/ FINALS/ A1 (20,16) , VI (20, 16) , A2 (20, 16) , V2 (20 , 16) , 

+  A3 (20,16) , V3 (20, 16) 

COMMON/ 10/  IN , I OUT , I DEBUG 
CHARACTER* 12  FNAME 
CHARACTER* 4  NAMEP 
CHARACTER* 1  YESNO 
REAL* 8  RNUM 

COMPLEX  CFAC , CFAC1 , CFAC2 , CFAC3 , CFAC4 , CF 
COMPLEX  CBLKF, CDEN, CH,CC, CM 
COMPLEX  CBLKR, CBLKB, CSHRB 

DOUBLE  PRECISION  TCQ1R,TCQ1C,TCQ2R,TCQ2C,TCQ3R,TCQ3C,TDSCRR 
DOUBLE  PRECISION  TDSCRC , TROOTR , TROOTC , TDEN , TNUMR , TNUMC , TLSQR 
DOUBLE  PRECISION  TLSQC , TLR , TLC , THOLD , DIV 

DIMENSION  RARR(100, 2) ,VARR(100,2) ,SARR(100,2) 

DIMENSION  DSARR ( 100 , 2 ) , DEARR( 100 , 2 ) , PARR (100, 2) 

A  RANDOM  NUMBER  GENERATOR  SEED  FOR  EACH 
OF  THE  SIX  VARIABLES  TO  BE  SAMPLED 

XI  -  .1 
X2  =  .2 
X3  =  .3 
X4  =>  .4 
X5  =  .  5 
X6  =  .  6 

PLUG-INS  FOR  THE  PRESENT 


INDICATE  PROVINCE  TO  BE  RUN 
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WRITE  (*,6000) 

FORMAT  (///'  ENTER  THE  FOUR  CHARACTER  NAME'/ 

•  OF  THE  PROVINCE  TO  BE  RUN. ') 

READ  (*,6001)  NAMEP 
FORMAT  (A4) 

WRITE  (*,6002)  NAMEP 

FORMAT  (/'  IS  ' , A4 , '  CORRECT  ?  (Y/N) » ) 

READ  (*,6003)  YESNO 
FORMAT  (Al) 

IF  (YESNO  .NE.  'Y*  .AND.  YESNO  .NE.  'y')  GOTO  100 


Ss 


M 


jl 
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READ  DATA  FROM  FILES 

WRITE  (FNAME, 6004)  NAMEP 
FORMAT  ( A4 , ' VOID . DAT • ) 

OPEN  (20,  FILE*FNAME , ACCESS= ' SEQUENTIAL ' , 

FORM= ' FORMATTED ' , STATUS* ' OLD ' ) 

READ  (20,6020)  NV, (VARR(J, 1) ,VARR(J, 2) , J=1,NV) 
FORMAT (IX, 13 , 100 (F7 . 4 , F7 . 3) ) 

CLOSE  (20,STATUS='KEEP' ) 


WRITE  (FNAME, 6005)  NAMEP 
FORMAT  ( A4 , ' SGGR . DAT ' ) 

OPEN  (20, FILE=FNAME , ACCESS* 1  SEQUENTIAL ' , 

FORM* ' FORMATTED • , STATUS* ' OLD  * ) 

READ  (20,6020)  NR, (RARR( J, 1) , RARR ( J, 2 ) , J=1 , NR) 
CLOSE  (20, STATUS* 1  KEEP ' ) 

WRITE  (FNAME, 6006)  NAMEP 
FORMAT  (A4 , ' SZRO . DAT ' ) 

OPEN  ( 20 , FILE=FNAME , ACCESS* ' SEQUENTIAL ' , 

FORM* ' FORMATTED ' , STATUS* ' OLD ' ) 

READ  (20,6021)  NS, (SARR(J, 1) ,SARR(J, 2) , J=1,NS) 
FORMAT  (1X,I3,100(F7.4,E11.4) ) 

CLOSE  (20, STATUS* ' KEEP 1 ) 


>.V, 
->v  , 
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WRITE  (FNAME, 6007)  NAMEP 
FORMAT  (A4, 'POIS.DAT' ) 

OPEN  ( 2 0,FILE*FNAME, ACCESS* 'SEQUENTIAL' , 

FORM* ' FORMATTED ' , STATUS* ' OLD ' ) 

READ  (20,6022)  NP, (PARR( J, 1) , PARR ( J , 2 ) , J=1 , NP) 
FORMAT  ( IX, I3,100(F7.4,F5.3) ) 

CLOSE  ( 20, STATUS*' KEEP ' ) 

WRITE  (FNAME, 6008)  NAMEP 
FORMAT  ( A4 , ' DECS . DAT ' ) 

OPEN  (20, FILE* FNAME, ACCESS*' SEQUENTIAL' , 

FORM* ' FORMATTED ' , STATUS* ' OLD ' ) 

READ  (20,6022)  NDS, (DSARR(J, 1) , DSARR(J, 2) , J=1,NDS) 
CLOSE  ( 20, STATUS*' KEEP ' ) 

WRITE  (FNAME, 6009)  NAMEP 
FORMAT  (A4 ,  'DECE.DAT') 

OPEN  (20 , FILE=FNAME , ACCESS* ' SEQUENTIAL ' , 

FORM* ' FORMATTED ' , STATUS* ' OLD ' ) 

READ  (20,6022)  NDE , (DEARR ( J , 1) , DEARR ( J , 2 ) , J=1 , NDE) 
CLOSE  ( 20, STATUS*' KEEP ' ) 
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WRITE  (FNAME, 6010)  NAMEP 

6010  FORMAT  (A4 BIOT. DAT ' ) 

OPEN  ( 10 , FILE=FNAME , ACCESS= ' SEQUENTIAL ' , 

&  FORM= ' UNFORMATTED • , STATUS= ' NEW ' ) 

WRITE  (FNAME, 6011)  NAMEP 

6011  FORMAT  (A4 , 1 BIOT.MTX ' ) 

OPEN  ( 11 , FILE=FNAME , ACCESS= ' SEQUENTIAL ' , 

&  FORM=’ FORMATTED' ,STATUS=' NEW' ) 

INPUT  THE  NUMBER  OF  TIMES  TO  STEP 
THROUGH  THE  MONTE  CARLO  LOOP. 

150  WRITE  (*,6030) 

5030  FORMAT  (///'  ENTER  THE  NUMBER  OF  MONTE  CARLO  STEPS.’) 
CALL  ISREAL  ( 10 , RNUM , IER) 

IF  (IER  .NE.  0)  GOTO  150 
MCLOOP  =  RNUM 


-  MAIN  PROGRAM  - 

REQUESTING  THE  MODEL  CONSTANTS 
CALL  MODCONSTS  (ISED) 

REQUESTING  TESTING  DEPTHS  (UP  TO  20  ALLOWED) 

CALL  QDEPTH  (.FALSE.) 

DO  1  I  =  2,20 

IF  ( EVERY ( I ) . LE . 0 . 0 )  GOTO  2 

1  CONTINUE 
1*1  +  1 

2  NDEPTH  *  I  -  1 

REQUESTING  TESTING  FREQUENCIES  (UP  TO  16  ALLOWED) 

CALL  QFREQ  (.FALSE.) 

DO  3  I  =  2,16 

IF  ( FREQ ( I ) . LE . 0 . 0 )  GOTO  4 

3  CONTINUE 
1  =  1  +  1 

4  NFREQ  =1-1 

BEGINNING  OF  THE  MONTE  CARLO  LOOPING 

MIDDEP  =  NDEPTH  /  2  +  1 

MIDFRE  =  NFREQ  /  2  +  1 

INT  =  FLOAT (MCLOOP  +  79)  /  FLOAT(80) 

MCEND  =  MCLOOP  /  INT 

DO  99  M  =  1, MCLOOP 


INITIALIZATION 

CFAC  =  CMPLX(0.0, 1.0) 

EPRESS  =0.0 

RANDOM  SAMPLING  OF  VALUES  FOR  RHOR, 
VOID,  SO,  DECBS,  DECBE ,  AND  PRATIO 


non  ooo  ooo  ono  ooo  o  o  o  o  o  o  o  o  ooo  ooo  ooo  ooo  oooo  ooo 


WRITE  (*,7001) 

7001  FORMAT  ('  CALLS  RNDSMP  (FIRST)') 
CALL  RNDSMP  (RARR,NR, X1,RH0R) 

WRITE  (*,7000) 

7000  FORMAT  ('  RETURN') 

WRITE  (*,7002) 

7002  FORMAT  ('  CALLS  RNDSMP  (SECOND)') 
CALL  RNDSMP  (VARR, NV, X2 , VOID) 

WRITE  (*,7000) 

WRITE  (*,7003) 

7003  FORMAT  ( '  CALLS  RNDSMP  (THIRD) ' ) 
CALL  RNDSMP  (SARR, NS , X3 , SO ) 

WRITE  (*,7000) 

WRITE  (*,7004) 

7004  FORMAT  (*  CALLS  RNDSMP  (FOURTH)') 
CALL  RNDSMP  (DSARR, NDS , X4 , DECBS) 

WRITE  (*,7000) 

WRITE  (*,7005) 

7005  FORMAT  ('  CALLS  RNDSMP  (FIFTH)') 
CALL  RNDSMP  (DEARR,NDE , X5 , DECBE) 

WRITE  (*,7000) 

WRITE  (*,7006) 

7006  FORMAT  ('  CALLS  RNDSMP  (SIXTH)') 
CALL  RNDSMP  (PARR, NP, X6 , PRATIO) 

WRITE  (*,7000) 


BETA  =  POROS  (VOID) 

PERM  =  PMBTY  ( BETA , H0VEMK , S 0 ) 

PSIZE  =  PORSIZ  (BETA, SO) 

DO  10  I  =  1 , NDEPTH 
WRITE  (*,7017)  I, NDEPTH 
7017  FORMAT  ('  I  =  ',12,'  NDEPTH  =  *,T2) 

WRITE  (*,7007) 

7007  FORMAT  ('  CALLS  OVERPRES ' ) 

CALL  OVERPRES  (I) 

WRITE  (*,7000) 

WRITE  (*,7008) 

7008  FORMAT  ('  CALLS  FSHEAR ' ) 

CALL  FSHEAR 

WRITE  (*,7000) 

WRITE  (*,7009) 

7009  FORMAT  ('  CALLS  FSHRIAM ' ) 

CALL  FSHRIAM 

WRITE  (*,7000) 

WRITE  (*,7010) 

7010  FORMAT  ('  CALLS  FBULK ' ) 

CALL  FBULK 

WRITE  (*,7000) 

WRITE  (*,7011) 

7011  FORMAT  ('  CALLS  FBLKIAM ' ) 

CALL  FBLKIAM 

WRITE  (*,7000) 

WRITE  (*,7012) 

7012  FORMAT  ('  CALLS  BIOTCOF') 

CALL  BIOTCOF 

C  WRITE  (*,7000) 
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DO  20  J  =  1 , NFREQ 
WRITE  (*,7018)  J, NFREQ 
7018  FORMAT  ('  J  =  ',12,'  NFREQ  =  ',12) 

WRITE  (*,7013) 

7013  FORMAT  ('  CALLS  FRQVARS ' ) 

CALL  FRQVARS  (J) 

WRITE  (*,7000) 

WRITE  (*,7014) 

7014  FORMAT  ('  CALLS  SKLVN 1 ) 

CALL  SKLVN 

WRITE  (*,7000) 

WRITE  (*,7015) 

7015  FORMAT  ('  CALLS  BDISREL' ) 

CALL  BDISREL  (I,J) 

WRITE  (*,7000) 

20  CONTINUE 
10  CONTINUE 

CALL  GRFDAT ( 1 1 , INT , MCEND , M , MI DDE P , MI DFRE ) 

CALL  WRIT  (10, M,NDEPTH, NFREQ) 

99  CONTINUE 

CLOSE  (10) 

CLOSE  (11) 

STOP 

END 


MCPHYS 


END  OF  MAIN  PROGRAM 


SUBROUTINES  USED  BY  THE  PROGRAM 


SUBROUTINE  MODCONSTS  (ISED) 


THIS  SUBROUTINE  PROCEEDS  TO  FILL  THE  VARIOUS  CONSTANTS 
USED  IN  THIS  VERSION  OF  BIOT-STOLL. 


COMMON/ CONSTS/ BLKRR , RHOF , BLKFR , ETA , ALPHA , HOVEMK , 
+  GC1 , GC2 , GC3 , VDH 

CHARACTER*!  ANS 


BLKRR=4 . 2E11 
RHOF®1.025 
BLKFR®0 . 2384E11 
ETA=0 . 018 
ALPHA® 1.2 5 
HOVEMK® 5 . 0 
GC2=2 . 97 
GC3=0 . 5 
ISED  =  10 
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98  WRITE  (*,100)  BLKRR 

100  FORMAT ( 1  Grain  Bulk  Modulus  =  ',E10.3,'  DYNES/CMA2 ' ) 

WRITE  (*,101)  RHOF 

101  FORMAT ('  Fluid  Density  =  ',F5.3,'  GRAMS/CMA3 • ) 

WRITE  (*,102)  BLKFR 

102  FORMAT ( '  Fluid  Bulk  modulus  =  ' ,E10.4,'  DYNES/ CM A 2  * ) 

WRITE  (*,103)  ETA 

103  FORMAT ('  Pore  Fluid  Viscosity  =  ' ,E10.3,'  POISE1) 

WRITE  (*,104)  ALPHA 

104  FORMAT ('  Structure  Factor  =  ' ,F5.2) 

WRITE  (*,105)  HOVEMK 

105  FORMAT ('  Kozeny-Carmen  Constant  *  ' ,F5.1) 

IF  (ISED  .EQ.  10)  THEN 

WRITE  (*,108) 

108  FORMAT  ('  Stoll  Stress  Formula  is  for  Silts/Clays. ' ) 

ELSE  IF  (ISED  .EQ.  1)  THEN 

WRITE  (*,109) 

109  FORMAT  ('  Stoll  Stress  Formula  is  for  Sand.') 

END  IF 

99  WRITE  (*,106) 

106  FORMAT ( '  Do  you  wish  to  change  any  or  all  of  these  ?  (Y/N):'/) 
READ  (*,107)  ANS 

107  FORMAT  (Al) 

IF  (ANS.EQ.'Y')  GOTO  1000 
IF  (ANS.EQ.'N')  GOTO  3000 
GOTO  99 
C 

1000  WRITE  (*,2001) 

WRITE  (*,1001) 

WRITE  (*,1002)  BLKRR 

1001  FORMAT ('  Give  a  Grain  Bulk  Modulus,  Real  Part', 

+  •  or  "RETURN"  to  continue.') 

1002  FORMAT ( '  (DEFAULT  =  ',E10.3,'  DYNES/CMA2) :  '/) 

READ  (* , 2001 , ERR=1000)  TEST1 

2001  FORMAT (E10 . 3 ) 

IF  (TEST1.NE.0.0)  BLKRR=TEST1 
C 

WRITE  (*,2002) 

1003  WRITE  (*,1004) 

WRITE  (*,1005)  RHOF 

1004  FORMAT ('  Give  a  Pore  Fluid  Density', 

+  '  or  "RETURN"  to  continue.') 

1005  FORMAT ( '  (DEFAULT  =  'F5.3,'  GRAMS/CMA3) :  '/) 

READ  (* , 2002 , ERR=1003 )  TEST2 

2002  FORMAT (F5. 3) 

IF  (TEST2.NE.0.0)  RHOF=TEST2 
C 

WRITE  (*,2003) 

1006  WRITE  (*,1007) 

WRITE  (*,1008)  BLKFR 

1007  FORMAT ('  Give  a  Fluid  Bulk  Modulus', 

+  '  or  "RETURN"  to  continue.') 

1008  FORMAT ('  (DEFAULT  =  ',E10.4,»  DYNES/CMA2) :  '/) 

READ  ( * , 2003 , ERR=1006)  TEST3 

2003  FORMAT (E10. 4) 

IF  (TEST3.NE.0.0)  BLKFR=TEST3 
C 

WRITE  (*,2004) 

1009  WRITE  (*,1010) 

WRITE  (*,1011)  ETA 
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1010  FORMAT ('  Give  a  Pore  Fluid  Viscosity', 

+  '  or  ' • RETURN' '  to  continue.') 

1011  FORMAT ('  (DEFAULT  =  ',E10.3,'  POISE):  '/) 

READ  ( * , 2004 ,ERR=1009)  TEST4 

2004  FORMAT (E10. 3) 

IF  (TEST4.NE.0.0)  ETA=TEST4 

C 

WRITE  (*,2005) 

1012  WRITE  (*,1013) 

WRITE  (*,1014)  ALPHA 

1013  FORMAT ('  Give  a  Structure  Factor', 

+  '  or  ''RETURN''  to  continue.') 

1014  FORMAT ( '  (DEFAULT  =  ',F5.2,'  }:  '/) 

READ  (* , 2005 , ERR=1012)  TEST5 

2005  FORMAT (F5. 2) 

IF  ( TEST5 . NE .0.0)  ALPHA=TEST5 

C 

WRITE  (*,2006) 

1015  WRITE  (*,1016) 

WRITE  (*,1017)  HOVEMK 

1016  FORMAT ('  Give  a  Kozeny-Carmen  Constant', 

+  '  or  ' ' RETURN ' '  to  continue . ' ) 

1017  FORMAT ( '  (DEFAULT  =  ',F5.1,'  ):  ’/) 

READ  (*, 2006 , ERR=1015)  TEST6 

2006  FORMAT (F5.1) 

IF  (TEST6.NE.0.0)  HOVEMK=TEST  6 
WRITE  (*,111) 

111  FORMAT  ('  Indicate  Stoll  Stress  Formula.'/ 

&  '  Default  is  Silts/Clays. '/ 

&  '  Do  you  want  to  change  it  to  Sand  Formula  ?  (Y/N) : ' ) 

READ  (*,107)  ANS 

IF  (ANS  .EQ.  'Y'  .OR.  ANS  .EQ.  'y')  ISED  =  1 
C 

GOTO  98 
C 

3000  IF  ( ISED . GE . 1 . AND . ISED . LE . 4 )  THEN 
GC1=1230.0 
VDH=0 . 5 
END  IF 
C 

IF  (ISED.GE. 5. AND. ISED. LE. 10)  THEN 
GC1=1630.0 
VDH=1 . 0 
END  IF 
C 
C 

WRITE  (*,2006) 

WRITE  (*,2006) 

WRITE  (*,2006) 

C 

RETURN 

END 

C 

C 

C 

SUBROUTINE  QDEPTH(LOOK) 

C 

C  ***  UPDATES  'EVERY'  FOR  DEPTHS 
C 

LOGICAL  LOOK 
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INTEGER  CASE, MAX 
CHARACTER  D*5,M*6 

COMMON  /VDEPTH/EVERY (20) 

D= ' DEPTH ' 

M-' METERS' 

•% 

EVERY  (1)  =0.0 
EVERY { 2 ) =0 . 5 
EVERY (3) =1. 0 
EVERY (4) =2.0 
EVERY (5) =5.0 
EVERY (6) =10.0 
EVERY ( 7 ) =2  0 . 0 
DO  1  1=8,20 

EVERY (I) =0 . 0 
1  CONTINUE 

■» 

5  IF  (.NOT. LOOK)  THEN 
CASE=MENU1 ( D , 5 ) 

ELSE 

CASE=1 
END  IF 
MAX=20 

IF  (CASE.EQ. 1)  THEN 

CALL  RLIST(D, 5,M, 6, MAX, EVERY) 

ELSE  IF  (CASE.EQ. 2)  THEN 

CALL  RADD(D, 5,M, 6, MAX, EVERY) 

ELSE  IF  (CASE.EQ. 3)  THEN 

CALL  RCHNGE (D, 5,M, 6, MAX, EVERY) 

ELSE  IF  (CASE.EQ. 4)  THEN 

CALL  RINSRT (D, 5 ,M, 6 , MAX, EVERY) 

ELSE  IF  (CASE.EQ. 5)  THEN 

CALL  RDELET ( D , 5 , M , 6 , MAX , EVERY ) 

END  IF 

IF  ( ( CASE . GT . 1  .AND.  CASE. LT. 6)  .OR. 

+  (CASE.EQ. 1  .AND.  .NOT. LOOK))  GOTO  5 
RETURN 
END 
C 

c 

SUBROUTINE  QFREQ  (LOOK) 

C 

C  ***  UPDATES  'FREQ'  FOR  FREQUENCIES 
C 

LOGICAL  LOOK 
INTEGER  CASE, MAX 
CHARACTER  F*9,H*5 
C 

COMMON  /VFREQ/FREQ ( 16 ) 

C 

F= ' FREQUENCY ' 

H= ' HERTZ ' 

FREQ (1) =100.0 
FREQ(2) =200.0 
FREQ (3) =400.0 
FREQ (4) =800.0 


C 
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FREQ (5) =1600.0 
DO  1  1=6,16 
FREQ (I) =0 
1  CONTINUE 
C 

5  IF  (.NOT. LOOK)  THEN 
CASE=MENU1 (F, 9) 

ELSE 

CASE=1 
END  IF 
MAX=16 

IF  (CASE.EQ. 1)  THEN 

CALL  RLIST ( F , 9 , H , 5 , MAX , FREQ ) 

ELSE  IF  (CASE.EQ. 2)  THEN 

CALL  RADD ( F , 9 , H , 5 , MAX , FREQ ) 

ELSE  IF  (CASE.EQ. 3)  THEN 

CALL  RCHNGE (F, 9 ,H, 5 ,MAX, FREQ) 

ELSE  IF  (CASE.EQ. 4)  THEN 

CALL  RINSRT ( F , 9 , H , 5 , MAX , FREQ ) 

ELSE  IF  (CASE.EQ. 5)  THEN 

CALL  RDELET ( F , 9 , H , 5 , MAX , FREQ ) 

END  IF 

IF  ( ( CASE . GT . 1  .AND.  CASE . LT . 6)  .OR. 

+  (CASE.EQ. 1  .AND.  .NOT. LOOK))  GO  TO  5 
RETURN 
END 

FUNCTION  QUERY (TOTAL) 

***  QUERY  PROMPTS  FOR  TOTAL  SELECTION 

INTEGER  TOTAL 

5  WRITE  (*,10)  TOTAL 
WRITE  (*,20) 

10  FORMAT  ('  ENTER  A  NUMBER  (1  -',12,')  CORRESPONDING  TO  THE') 
>  FORMAT  ('  OPTION  DESIRED.  (ENTER  0  FOR  MENU):  '/) 

READ  (* , 30 , ERR=5)  NUMBER 
30  FORMAT (12) 

IF  ( NUMBER. LT.O. OR. NUMBER. GT. TOTAL)  GOTO  5 
QUERY =NUMBER 
RETURN 
END 


FUNCTION  MENU1 (OPTION, LENGTH) 

***  MENU1  PRINTS  THE  DEPTHS  AND  FREQUENCY  MENUS 

CHARACTER  OPTION* 9 

IF  (LENGTH. EQ. 5)  THEN 
WRITE ( * , 3 )  OPTION 

3  FORMAT  ('O' , A5) 

ELSE 

WRITE ( * , 4 )  OPTION 

4  FORMAT ( *  0  * , A9 ) 
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END  IF 

5  CHOICE-QUERY (6) 

IF  (CHOICE. EQ.O)  THEN 
WRITE (*, 10) 


10 


+ 

+ 

+ 

+ 

+ 

+ 


FORMAT  ( 


OOPTIONS : 
1  — 
2  — 

3  — 

4  — 

5  — 

6  — 


GO  TO  5 


ELSE 


MENU 1-CHOICE 
END  IF 
RETURN 
END 


V 

LIST'/ 

ADD'/ 

CHANGE'/ 

INSERT'/ 

DELETE'/ 

EXIT  ' ) 


SUBROUTINE  RLIST ( PARAM , LP , UNITS , LU , MAX , ARRAY ) 


***  LISTS  REAL  'ARRAY'. 


CHARACTER  PARAM*9 , UNITS  *  6 
DIMENSION  ARRAY (MAX) 

COMMON  /VDEPTH/EVERY(20) 

COMMON  /VFREQ  /FREQ (16) 

DO  5  I— 1 , MAX 

IF  (MAX.EQ.20)  THEN 
ARRAY (I) -EVERY (I) 

ELSE 

ARRAY (I) -FREQ (I) 

END  IF 
5  CONTINUE 

IF  (MAX.EQ.20)  THEN 

WRITE (*,30)  1, ARRA( 1) 

ELSE 

WRITE (* , 35)  1 , ARRAY ( 1) 

END  IF 

30  FORMAT  ('  ',12,':  \F6.2) 

35  FORMAT  ('  ',12,':  ',G10.5) 

DO  40  1-2, MAX 

IF  (ARRAY (I) . NE . 0 . )  THEN 
IF  (MAX.EQ.20)  THEN 

WRITE (*,30)  I, ARRAY (I) 

ELSE 

WRITE (*, 35)  I , ARRAY ( I ) 
END  IF 

ELSE 

GO  TO  41 
END  IF 

40  CONTINUE 

41  CONTINUE 
RETURN 
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SUBROUTINE  RADD ( PARAM , LP , UNITS , LU , MAX , ARRAY ) 


C 

C  *** 

C 


***  ADD  APPENDS  DATA  AFTER  THE  LAST  VALID  ELEMENT  OF  'ARRAY' 

COMMON  / VDEPTH/EVERY (20) 

COMMON  /VFREQ  /FREQ (16) 

CHARACTER  PARAM*9,UNITS*6 
DIMENSION  ARRAY (MAX) 

IF  (ARRAY (1) .EQ.O.  .AND.  ARRAY (2) . EQ. 0 . )  THEN 
NDX-1 

ELSE 

DO  10  1-2, MAX 

IF  (ARRAY(I) .EQ.O.)  THEN 
NDX-I-1 
GO  TO  11 
END  IF 

10  CONTINUE 

11  CONTINUE 
END  IF 

12  IF  (MAX.EQ.20)  THEN 

WRITE (*,20)  ARRAY (NDX) 

ELSE 

WRITE (*,25)  ARRAY (NDX) 

END  IF 

15  WRITE  (*,30) 

20  FORMAT  ('OLAST  IS:  ',F6.2) 

25  FORMAT  ('OLAST  IS:  ',G8.3) 

30  FORMAT  ( ' $ENTER  ADDITION (S) ;  (0  TO  EXIT)  ') 

40  WRITE (*, 42) 

42  FORMAT  ('  :  '/) 


IF  (MAX.EQ.20)  THEN 

■i 

READ(*,45,ERR-12)  VALUE 

t  £ 

ELSE 

READ(* , 47 , ERR—12)  VALUE 

ENDIF 

1  V 

45 

FORMAT  (F6.2) 

1  >. 

47 

FORMAT  (G8.0) 

i  v 

IF  (VALUE. LT.O.)  THEN 

.. 

WRITE (*,50) 

50 

FORMAT  ('  ENTER  A  POSITIVE 

t  V 
« 

GO  TO  40 

1* 

1 

ELSE  IF  (VALUE . NE . 0 . )  THEN 

» 

NDX— NDX+1 

<•  .** 

: 

ARRAY (NDX) -VALUE 

w 

GO  TO  40 

END  IF 

S  ry, 

1  *  m 

DO  60  I— 1 ,MAX 

■ 

IF  (MAX.EQ.20)  THEN 

V 

; 

EVERY (I) -ARRAY (I) 

ELSE 

:•  -v 

FREQ (I) -ARRAY (I) 

■, 

END  IF 

60 

CONTINUE 

i 

RETURN 

SUBROUTINE  RCHNGE ( PARAM , LP , UNITS , LU , MAX , ARRAY ) 
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CHANGES  VALUE  IN  SPECIFIED  INDEX  OF  REAL  'ARRAY1 

CHARACTER  PARAM*9 , UNITS* 6 , ANSI* 1 , ANS2*1 
DIMENSION  ARRAY (MAX) 

COMMON  /VDEPTH/EVERY (20) 

COMMON  /VFREQ  /FREQ (16) 


WRITE (* , 5) 

FORMAT ( ' 0 ' ) 

WRITE (*, 20) 

FORMAT  ('  INDEX  OF  ONE  TO  CHANGE;  (0  TO  EXIT): 
READ (*,35, ERR=10 )  NDX 
FORMAT  (12) 

IF  (NDX.EQ.O)  GO  TO  1000 
IF  (NDX. GT. MAX  .OR.  NDX.LT.O)  THEN 
WRITE (*,40)  MAX 

FORMAT  ('  INDEX  MUST  BE  A  NUMBER,  1- * , 12 , ' 
GO  TO  10 
END  IF 

IF  (MAX. EQ. 20)  THEN 

WRITE (* , 50)  PARAM, NDX, ARRAY (NDX) , UNITS 

ELSE 

WRITE ( * , 52 )  PARAM, NDX, ARRAY (NDX) , UNITS 
END  IF 

FORMAT  ('  ' ,A6, '  #  ',12,*  IS  ',F6.2,'  ',^6,'.') 
FORMAT  ('  ' , A9 , 1  #  ',12,'  IS  *,G8.3,'  ',A6,'.') 
WRITE (*,60) 

FORMAT  ('  DO  YOU  WANT  TO  CHANGE  IT?  (Y/N) :  '/) 
READ (*,998)  ANSI 
IF  (ANS1.EQ. • Y * )  THEN 
WRITE (*,80) 

FORMAT  ('  ENTER  NEW  VALUE:  '/) 

IF  (MAX.EQ.20)  THEN 

READ (*,90, ERR=7  0 )  VALUE 

ELSE 

READ (*,95, ERR=7  0 )  VALUE 
END  IF 

FORMAT  (F6.2) 

FORMAT  (G8.0) 

IF  (VALUE.LT. 0.)  THEN 
WRITE (*, 100) 

FORMAT  ( ' $ ENTER  A  NON-NEGATIVE  NUMBER 
GO  TO  70 
END  IF 

ARRAY (NDX) =VALUE 
ELSE  IF  (ANS1.NE. ' N  * )  THEN 
GO  TO  55 
END  IF 
WRITE (*, 120) 

FORMAT  ('  CHANGE  ANOTHER  ?  (Y/N):  '/) 

READ (*,998)  ANS2 
IF  (ANS2.EQ. 'Y' )  THEN 
GO  TO  10 

ELSE  IF  (ANS2 .NE. 'N' )  THEN 
GO  TO  110 
END  IF 

DO  130  1=1, MAX 

IF  (MAX.EQ.20)  THEN 
EVERY ( I ) = ARRAY ( I ) 
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ELSE 

FREQ ( I ) = ARRAY ( I ) 
END  IF 


C  ***  INSERTS  VALUES  BEFORE  THE  SPECIFIED  INDEX  IN  REAL  'ARRAY' 
C 

CHARACTER  PARAM*9 ,UNITS*6 , ANS1*1 
DIMENSION  ARRAY (MAX) 

LOGICAL  L 

COMMON  /VDEPTH/EVERY (20) 

COMMON  /VFREQ  /FREQ (16) 


WRITE ( * ,  5 ) 

5  FORMAT  ('0') 

10  WRITE (*,20) 

20  FORMAT  ('  INDEX  OF  VALUE  TO  INSERT  BEFORE;  (0  TO  EXIT) 
READ(*,  25,ERRS*10)  NDX 
25  FORMAT  (12) 

L-(NDX.EQ.O  .OR.  (NDX.EQ.l  .AND.  ARRAY ( 1) .EQ. 0 . 0) ) 

IF  (L)  GO  TO  1000 

IF  (NDX. GE. MAX  .OR.  NDX.LT.O)  THEN 
WRITE (*,30)  MAX-1 

30  FORMAT  ('  INDEX  MUST  BE  A  NUMBER,  1-',I2,';  '/) 

GO  TO  10 
END  IF 

40  IF  (MAX. EQ. 20)  THEN 

WRITE  (*,50)  NDX, ARRAYNDX) , UNITS 

ELSE 

WRITE  (*,55)  NDX, ARRAY (NDX) , UNITS 
END  IF 

50  FORMAT  ( '  DO  YOU  WANT  TO  ENTER  VALUES  BEFORE  #  '  , 12 , 


'/) 


' , F6 . 2 , '  1 , A6, ' ?  (Y/N):  •/) 


55  FORMAT  ('  DO  YOU  WANT  TO  ENTER  VALUES  BEFORE  #',I2, 


60 

70 

80 

84 


READ (*,998)  ANSI 
IF  (ANS1.EQ. ' Y' )  THEN 


'  ,  G8 . 3  ,  '  \A6  ,*?  (Y/N):  •/) 


90 


95 

97 


WRITE (*,70) 

FORMAT ( '  ENTER  VALUE(S)  TO  INSERT,  (0  TO  EXIT) 
WRITE (*,84) 

FORMAT  ('  :  '/) 

IF  (NDX. GE. MAX)  THEN 
WRITE (*,90)  MAX 

FORMAT  ('  ONLY  ',12,'  VALUES  ALLOWED.') 

ELSE 

IF  (MAX. EQ. 20)  THEN 

READ ( * , 95 , ERR=60)  VALUE 

ELSE 

READ (* , 97 , ERR=60)  VALUE 
END  IF 

FORMAT  (F6.2) 

FORMAT  (G8.0) 
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130 

CONTINUE 

998 

FORMAT  ( 1A1) 

1000 

CONTINUE 

■■■V.y. 

RETURN 

£  -  i 

END 
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SUBROUTINE  RINSRT ( PARAM , LP , UNITS , LU , MAX , ARRAY ) 

V. <■, 

K.-; 


I _ 


•  V  "«* 


&8! 

i 


-■v-> 


l 

•  •  • 


.vvi 


A 


.  /.  -•  •• 


\ 

V 


.-V-V-V. 


c 

c 

c 


IF  (NDX.GT.l  .AND.  VALUE . EQ . 0 . 0 )  GO  TO  1000 

* 

*  FIND  INDEX  OF  LAST  ARRAY  VALUE 

* 

DO  100  1=2, MAX 

IF  (ARRAY ( I ) . EQ . 0 . )  THEN 
LAST=I 
GO  TO  101 
END  IF 

100  CONTINUE 

101  CONTINUE 

C  * 

C  *  SHIFT  RIGHT  &  INSERT  VALUE 

C  * 

DO  110  I=LAST,NDX,-1 

ARRAY ( 1+ 1 ) =ARRAY ( I ) 

110  CONTINUE 

ARRAY (NDX) =VALUE 
NDX=NDX+1 
GO  TO  80 
END  IF 

ELSE  IF  (ANS1.NE. * N * )  THEN 
GO  TO  40 
END  IF 

DO  120  1=1, MAX 

IF  (MAX.EQ.20)  THEN 
EVERY ( I ) =ARRAY ( I ) 

ELSE 

FREQ ( I ) =ARRAY ( I ) 

END  IF 
120  CONTINUE 
998  FORMAT  (1A1) 

1000  CONTINUE 
RETURN 
END 
C 

C 

SUBROUTINE  RDELET ( PARAM , LP , UNITS , LU , MAX , ARRAY ) 

C 

C  ***  DELETES  RANGE  OF  VALUES  BETWEEN  SPECIFIED  INDECES 
C 

CHARACTER  PARAM*9 ,UNITS*6 , ANS1*1 , ANS2*1 
DIMENSION  ARRAY (MAX) 

COMMON  /VDEPTH/EVERY (20) 

COMMON  /VFREQ  /FREQ (16) 

C 

WRITE ( * , 5 ) 

5  FORMAT  ('0') 

10  WRITE ( * , 20) 

20  FORMAT  ('  FIRST  INDEX  OF  RANGE  TO  DELETE;  (0  TO  EXIT) 
READ ( * , 25 , ERR=10)  NDX1 
25  FORMAT  (12) 

IF  (NDX1.EQ.0)  GO  TO  1000 
IF  (NDX1.GT.MAX)  GO  TO  10 
30  WRITE ( * , 40) 

40  FORMAT  ( ' $LAST  INDEX  OF  RANGE  TO  DELETE,  (0  TO  EXIT): 
READ (*,25, ERR=3  0 )  NDX2 
IF  (NDX2.EQ.0)  GO  TO  1000 

IF  (NDX2.LT. NDX1  .OR.  NDX2.GT.MAX)  GO  TO  30 


MCPHYS 


DO  50  I=NDX1,NDX2 

IF  (MAX.EQ.20)  THEN 

WRITE ( * , 60)  I , ARRAY ( I ) 

ELSE 

WRITE (* , 62)  I , ARRAY ( I ) 
END  IF 


(I 

50 

CONTINUE 

60 

FORMAT  ( •  • ,12, • :  • ,F6.2) 

62 

FORMAT  ( •  •  ,12, • :  • ,G8.3) 

,v 

65 

WRITE (*,70) 

S" 

70 

FORMAT  (•  DELETE?  (Y/N):  •/) 

READ(*, 998)  ANSI 

N=MAX-NDX2 

* 

IF  (ANS1.EQ. 'Y* )  THEN 

IF  (NDX1 . LT.MAX)  THEN 

DO  80  1=1, N 

ARRAY (-1+NDX1+I) : 

80 

CONTINUE 

IF  (NDX1.EQ.NDX2)  N=] 
DO  85  I=NDX1+N , MAX 

ARRAY (I) =0.0 

.*•* 

85 

CONTINUE 

END  IF 

r  *. 

ARRAY (MAX) =0 . 0 

f:: 

ELSE  IF  (ANS1.NE. 'N' )  THEN 

GO  TO  65 

END  IF 

V 

DO  87  1=1, MAX 

998 

1000 


IF  (MAX.EQ.20)  THEN 
EVERY ( I ) = ARRAY ( I ) 

ELSE 

FREQ ( I ) =ARRAY ( I ) 

END  IF 
CONTINUE 
WRITE(*, 100) 

ORMAT  ('  DELETE  ANOTHER  VALUE  ?  (Y/N) :  '/) 
READ (*,998)  ANS2 
IF  (ANS2.EQ. 'Y' )  THEN 
GO  TO  10 

ELSE  IF  (ANS2.NE. 'N')  THEN 
GO  TO  90 
END  IF 
FORMAT  ( 1A1) 

CONTINUE 

RETURN 

END 


FUNCTION  POROS  (VOID) 

FUNCTION  TO  CALCULATE  POROSITY 
POROS  =  VOID/ ( 1+VOID) 


RETURN 

END 
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FUNCTION  PMBTY  (BETA, HOVEMK, SO) 

FUNCTION  TO  CALCULATE  PERMEABILITY 
USING  KOZENY -CARMEN 

PMBTY  =  BETA**3/ (HOVEMK*SO**2* (1-BETA) **2) 

RETURN 

END 


FUNCTION  PORSIZ  (BETA, SO) 

FUNCTION  TO  CALCULATE  PORE  SIZE 
USING  HOVEM  AND  INGRAM  FOR  PORE 
SIZE  BEING  EQUAL  TO  TWICE  THE 
HYDRAULIC  RADIUS  WHICH  IN  TURN 
DEPENDS  UPON  MEAN  SPECIFIC  SURFACE  (SO) 

PORSIZ  =  2*BETA/( (1-BETA) *S0) 

RETURN 

END 


SUBROUTINE  OVERPRES  (ICOUNT) 

SUBROUTINE  TO  CALCULATE  OVER  BURDEN  PRESSURE 

COMMON/ CONSTS/BLKRR , RHOF , BLKFR, ETA , ALPHA , HOVEMK , 
+  GC1 , GC2 , GC3 , VDH 

COMMON/VDEPTH/EVERY (20) 

COMMON/ DISTVALS/RHOR, VOID, SO , DECBS , DECBE , PRATIO 
COMMON/ FCTVLS/ BETA, PERM, F^TZE 
COMMON/DEPLOOP/EPRESS , G, C  T,E,EPRI 

G  =  GRAVITATIONAL  CON.  IT  (CM/ SEC) 

G  =  980.665 
IF  (ICOUNT. EQ. 1)  THEN 
DDEPTH  =  EVERY (1) 

ELSE 

DDEPTH  =  EVERY (ICOUNT)  -  EVERY ( I COUNT- 1) 

END  IF 

DPRESS  =  DDEPTH* (RHOR  -  RHOF)*(l  -  BETA) *G*100 
EPRESS  =  EPRESS  +  DPRESS 

RETURN 

END 


SUBROUTINE  FSHEAR 

SUBROUTINE  TO  CALCULATE  FRAME  SHEAR  MODULUS 
USING  THE  STOLL  STRESS  PROCEDURE  WHICH  INCLUDES 
A  DEPTH  DEPENDANCE 
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COMMON/ CONSTS/BLKRR , RHOF , BLKFR, ETA , ALPHA , HOVEMK , 
+  GC1 , GC2 , GC3 , VDH 

COMMON/DISTVALS/RHOR, VOID, SO, DECBS , DECBE, PRATIO 
COMMON/DEPLOOP/EPRESS , G , GPRI , E , EPRI 
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SGMVRT  -  EPRESS* (1.45E-05) 

SGMBAR  -  ((1.0  +  ( 2 . 0*VDH) ) *SGMVRT) /3 . 0 
TLC  -  GC2  -  VOID 

G  «  (GC1* (TLC*TLC) * (SGMBAR**GC3 ) ) / ( 1 . 0  +  VOID) 
G  -  G/ ( 1 . 45E-05) 


RETURN 

END 


SUBROUTINE  FSHRIAM 

SUBROUTINE  TO  CALCULATE  THE  IMAGINARY  PART  OF 
THE  FRAME  SHEAR  MODULUS 

COMMON/ DISTVALS/RHOR , VOID , SO , DECBS , DECBE , PRATIO 
COMMON/DEPLOOP/EPRESS , G , GPRI , E, EPRI 


PI  =  3.1415926535898 
GPRI  =  (G/PI) * DECBS 


RETURN 

END 


SUBROUTINE  FBULK 

SUBROUTINE  TO  CALCULATE  THE  FRAME  BULK  MODULUS 
AS  IT  IS  DERIVED  FROM  THE  FRAME  SHEAR  MODULUS 

COMMON/DISTVALS/RHOR, VOID, SO, DECBS , DECBE, PRATIO 
COMMON/DEPLOOP/EPRES , G , GPRI , E , EPRI 

E  «  ( (2.0*G) *(1.0  +  PRATIO) )/ (3.0* (1.0  -  (2 . 0*PRATIO) ) ) 

RETURN 

END 


SUBROUTINE  FBLKIAM 

SUBROUTINE  TO  CALCULATE  THE  IMAGINARY  PART  OF 
THE  FRAME  BULK  MODULUS 

COMMON/DISTVALS/RHOR, VOID , SO , DECBS , DECBE , PRATIO 
COMMON/DEPLOOP/EPRESS , G , GPRI , E , EPRI 

PI  -  3.1415926535898 
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EPRI  =  (E/PI) *DECBE 

RETURN 

END 


SUBROUTINE  BIOTCOF 

SUBROUTINE  TO  CALCULATE  THE  BIOT  COEFFICIENTS 
COMPLEX  CBLKR , CBLKF , CBLKB , CSHRB , CDEN , CFAC1 , CH , CC , CM 


i 

•c* 

v, 
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COMMON/ CONSTS/BLKRR ,  RHOF ,  BLKFR ,  ETA ,  ALPHA ,  HOVEMK ,  ..w 

+  GC1 ,  GC2 ,  GC3 ,  VDH  V 

COMMON/  FCTVLS/ BETA ,  PERM ,  PS  I Z  E 

COMMON/DEPLOOP/EPRESS  ;  G ,  GPRI ,  E ,  EPRI  . 

COMMON/ BCOEF/ CBLKR,  CBLKF ,  CBLKB ,  CSHRB ,  CDEN ,  CFAC1 ,  CH ,  CC ,  CM  v 


CBLKR  =  CMPLX ( BLKRR ,0.0) 
CBLKF  =  CMPLX ( BLKFR , 0 . 0 ) 
CBLKB  =  CMPLX (E, EPRI) 
CSHRB  =  CMPLX (G, GPRI) 


CDEN  =  (CBLKR* (1.0  +  (BETA* ( (CBLKR/CBLKF)  -  1.0))))  -  CBLKB 
CFAC1  =  CBLKR  -  CBLKB 

CH  =  ( (CFAC1*CFAC1)/CDEN)  +  CBLKB  +  ( (4 . 0/3 . 0) *CSHRB) 

CC  -  (CBLKR*CFAC1) /CDEN 
CM  =  (CBLKR* CBLKR) /CDEN 


RETURN 

END 


SUBROUTINE  FRQVARS  (ICOUNT) 
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SUBROUTINE  TO  CALCULATE  OMEGA,  EKAPPA,  AND  EM 


CC 

CC 


COMMON/ CONSTS/ BLKRR , RHO , BLKFR , ETA , ALPHA , HOVEMK , 
+  GC1 , GC2 , GC3 , VDH 

COMMON/ VFREQ/ FREQ ( 16 ) 

COMMON/ FCTVLS/ BETA, PERM, PSIZE 
COMMON/ FVARS/OMEGA , EKAPPA , EM 

PI  =  3.1415926535898 

OMEGA  =  ( 2. 0*PI ) *FREQ( ICOUNT) 

EKAPPA  =  PSIZE*SQRT( (OMEGA*RHOF) /ETA) 

EM  =  (ALPHA*RHOF) /BETA 

RETURN 

END 


V/ 


SUBROUTINE  SKLVN 

COMPUTES  CORRECTION  FACTOR  FOR  DYNAMIC  VISCOSITY  v 

SEE  ABRAMOWITZ  AND  STEGUN,  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS,  P. 
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COMMON/ BCOEF/CBLKR, CBLKF , CBLKB, CSHRB, CDEN , CFAC1 , CH, CC, CM  MCPHYS 

COMMON/ FVARS/OMEGA , EKAPPA , EM 

COMMON/ COMKLV/ CFAC , CFAC2 , CFAC3 , CFAC4 , CF 

COMPLEX  CFAC1 , CFAC2 , CFAC3 , CFAC4 , CFAC, CF, CPHIM, CPHIP, CTHTM , CTHTP 
COMPLEX  CBLKR , CBLKF , CBLKB , CSHRB , CDEN , CH , CC , CM , CTEMPA , CTEMPB 
C 

PI  «  3.1415926535898 
C 

IF (EKAPPA. GE. 119)  GO  TO  2003 
2000  IF ( EKAPPA . GE .8.0)  GO  TO  2005 
X=EKAPPA/8 . 0 
XSQ=X*X 
XFO=*XSQ*XSQ 

CFAC1=1 . 0  +  XFO* (-64 . 0  +  XFO* ( 113 . 77777774  +  XFO* (-32 . 36345652  + 

1  XFO* (2. 64191397  +  XFO* (-8 . 349609E-02  +  XFO* ( 1 . 22552E-03  + 

2  XFO* (-9 . 01E-06) )))))) 

CFAC2=XSQ* (16.0  +  XFO* (-113 . 77777774  +  XFO* (72 . 81777742  4- 

1  XFO* (-10. 56765779  +  XFO* ( 0 . 52185615  +  XFO* (-1 . 103667E-02  + 

2  XFO* (1. 1346E-04) )))))) 

CFAC3=(EKAPPA*XSQ) *(-4.00  +  XFO* (14 . 22222222  +  XFO* (-6 . 06814810  + 

1  XFO*(0. 66047849  +  XFO* (-2 . 609253E-02  +  XFO* (4 . 5957E-04  + 

2  XFO* (-3 . 94E-06) )))))) 

CFAC4=EKAPPA* (0.5  +  XFO* (-10 . 66666666  +  XFO* ( 11 . 37777772  + 

1  XFO*(-2. 31167514  +  XFO* (0. 14677204  +  XFO* (-3 . 79386E-03  + 

2  XFO* (4 . 609E-05) )))))) 

CF= (CFAC3  +  ( CFAC*CFAC4 ) ) / ( CFAC1  +  (CFAC*CFAC2) ) 

GO  TO  2010 

USE  ASYMPTOTIC  APPROXIMATION  FOR  LARGE  KAPPA 
KELVIN (KAPPA) =KAPPA/ ( 4  *SQRT ( 2 ) ) 

=.176777*KAPPA 
MARCH, 1984  DAW 

2003  CF=.176777*EKAPPA 
GO  TO  2099 


2005  X=8 . 0/ EKAPPA 
XSQ=X*X 

CFAC1=(0. 0,-0. 3926991)  +  XSQ* ( ( 0 . 0 , -9 . 765E-04 )  + 

1  XSQ* ( (-2 . 52E-05 ,0.0)  +  XSQ* ( (6 . 0E-07 , 1 . 9E-06) ) ) ) 
CFAC2=X*((1.10486E-02,-1.10485E-02)  +  XSQ* ( (-9 . 06E-05 , -9 . 01E-05)  + 
1  XSQ* ( (-3 . 4E-06 , 5 . IE-06) ) ) ) 

CTHTP=CFAC1  +  CFAC2 
CTHTM=CFAC 1  -  CFAC 2 

CFAC3=(0. 7071068, 0.7071068)  +  XSQ* ( ( -1 . 3813E-03 , 1 . 38 11E-03 )  + 

1  XSQ* ( (3.46E-05,3.38E-05)  +  XSQ* ( ( 1 . 6E-06 , -3 . 2E-06) ) ) ) 

CFAC4=»X*  (  (-6 . 25001E-02  ,  -1 . 0E-07 )  +  XSQ*  (  ( 5 . 0E-07 , 2 . 452E-04 )  + 

1  XSQ* ( ( 1 . 17E-05 , -2.4E-06) )  )  ) 

CPHIP=CFAC3  +  CFAC4 
CPHIM=CFAC3  -  CFAC4 
CFAC1=1 . 0/ (SQRT ( 2 . 0*PI*EKAPPA) ) 

CFAC2= (1.0  +  CFAC) / (SQRT (2.0) ) 

CFAC3= ( -CFAC2  *EKAPPA)  +  CTHTM 
CFAC4* ( CFAC2  *EKAPPA)  +  CTHTP 
X-4.0 

C  write  (*,*)  cexp (cfac3/x) , cfac3 , x 

CF*CEXP (CFAC3/X) 

C  write  (*,*)  cf 

CF-CF*CF 


c 


write  (*,*)  cf 
CF=CF*CF 
C  write  (*,*)  pi,cfacl,cf 

CFAC3= ( PI *CFAC1 ) *CF 

C  write  (*,*)  cexp(cfac4/x) ,cfac4,x 

CF=CEXP (CFAC4/X) 

C  write  (*,*)  cf 

CF=CF*CF 

C  write  (*,*)  cf 

CF=CF*CF 

C  write  (*,*)  cfacl,cf 

CFAC4=CFAC1*CF 
C  write  (*,*)  cfac,pi 

CFAC1=CFAC/PI 

C  write  (*,*)  '  ONE ' , cf ac4 , cphip , cf acl , cf ac3 , cphim 

CTEMPA  =  (CFAC4  *  CPHIP)  -  (CFAC1  *  CFAC3  *  CPHIM) 

C  write  (*,*)  '  A' , ctempa , cfac4 , cf acl , cfac3 

CTEMPB  =  CFAC4  +  (CFAC1  *  CFAC3 ) 

C  write  (*,*)  '  B',ctempb 

TEMPRA  =  REAL (CTEMPA) 

TEMPRB  =  REAL (CTEMPB) 

TEMPI B  =  AIMAG (CTEMPB) 

C  write  (*,*)  '  C',tempra 

TEMPI A  =  AIMAG (CTEMPA) 

TEMPRB  =  TEMPRB  /  TEMPRA 
TEMPI B  =  TEMPI B  /  TEMPRA 
TEMPIA  =  TEMPIA  /  TEMPRA 
TEMPRA  =  TEMPRA  /  TEMPRA 
CTEMPA  =  CMPLX (TEMPRA, TEMPIA) 

C  CTEMPA  =  CTEMPA  /  TEMPRA 

C  write  (*,*)  '  D', ctempa 

CTEMPB  =  CMPLX (TEMPRB, TEMPI B) 

C  CTEMPB  -  CTEMPB  /  TEMPRA 

C  write  ( * , * )  1  E ' , ctempb 
CF  =  CTEMPA  /  CTEMPB 

C  CF=( (CFAC4*CPHIP) - (CFAC1*CFAC3*CPHIM) ) / (CFAC4+ (CFAC1*CFAC3 ) ) 

C2010  write  (*,*)  '  TWO' ,ekappa,cf ,cfac 

2010  CF= (0.25  * (EKAPPA*CF) ) / ( 1 . 0  -  ((2.0  *CF) / (CFAC*EKAPPA) ) ) 

2099  RETURN 
END 
C 
C 
C 

SUBROUTINE  BDISREL  (I,J) 

CC 

CC  SOLVE  BIOT  DISPERSION  RELATIONS 

CC 

COMMON/ CONSTS/BLKRR , RHOF , BLKFR , ETA , ALPHA , HOVEMK , 

+  GC1 , GC2 , GC3 , VDH 

COMMON/DISTVALS/RHOR, VOID, SO, DECBS , DECBE, PRATIO 
COMMON/ FCTVLS/ BETA , PERM , PS I Z  E 

COMMON/BCOEF/CBLKR , CBLKF , CBLKB , CSHRB , CDEN , CFAC1 , CH , CC , CM 
COMMON/ FVARS/OMEGA , EKAPPA , EM 
COMMON/ COMKLV/CFAC, CFAC2 , CFAC3 , CFAC4 , CF 
COMMON/FINALS/A1 (20 , 16) ,V1(20,16) ,A2(20,16) ,V2(20,16) , 

+  A3 (20, 16) ,V3 (20, 16) 

CC 

COMPLEX  CFAC1 , CFAC2 , CFAC3 , CFAC4 , CFAC , CF 
COMPLEX  CBLKF , CDEN , CH , CC , CM , CQ1 , CQ2 , CQ3 
COMPLEX  CBLKR , CBLKB , CSHRB , CTEMPA , CTEMPB 


MCPHYS 


DOUBLE  PRECISION  TCQ1R,TCQ1C,TCQ2R,TCQ2C,TCQ3R,TCQ3C,TDSCRR 
DOUBLE  PRECISION  TDSCRC , TROOTR , TROOTC , TDEN , TNUMR , TNUMC , TLSQR 
DOUBLE  PRECISION  TLSQC , TLR , TLC , THOLD , DIV 


DIMENSION  THOLD (6) 


RHO  =  (RHOF*BETA)  +  (RH0R*(1.0  -  BETA)) 

NOFRAM  =  1 

D  =  10.0 

EX  =  CABS ( CBLKB ) 

GX  *  CABS (CSHRB) 

IF  (EX.LT.D.AND.GX.LT.D)  NOFRAM  =  2 

CC 

CC***  CALCULATE  VELOCITY  AND  ATTENUATION 
OSQ=OMEGA*OMEGA 
GO  TO  (330,320),  NOFRAM 
320  CFAC1= (CFAC*ETA*CF) / (OMEGA* PERM) 

CFAC2=( ( (ALPHA/BETA)  -  (RHOF/RHO)  )  *RHOF)  -  CFAC1 
CFAC3=RHO  +  (( (ALPHA/ BETA)  -  2.0)*RHOF)  -  CFAC1 
CFAC4=1.0/( (BETA/CBLKF)  +  ((1.0  -  BETA) /CBLKR) ) 

CF=( (RHO/CFAC4) * (CFAC2/CFAC3 ) ) *OSQ 
TLSQR=REAL(CF) 

TLSQC=AIMAG ( CF) 

CALL  SCROOT  (TLSQR, TLSQC, TLR, TLC) 

A1 (I , J) =-434 . 294481*TLC*2 . 0 
VI ( I , J) - (OMEGA/TLR) /100 . 0 
A2 (I , J) =0 . 0 
V2 (I , J) =0 . 0 
A3 (I , J) =0 . 0 
3 (I, J)=0. 0 
GO  TO  350 

330  CFACl=RHOF*OSQ 

CFAC2= (CFAC*OMEGA*CF*ETA) /PERM 

CFAC3=EM*OSQ 

CFAC4=RHO*OSQ 

CQ1=(CC*CC)  -  (CH*CM) 

CQ2=( (CH*CFAC3)  +  (CM*CFAC4))  -  (2 . 0*CC*CFAC1)  -  (CH*CFAC2) 
CQ3=(CFAC1*CFAC1)  -  (CFAC3*CFAC4)  +  (CFAC4*CFAC2 ) 

TCQ 1R=REAL ( CQ 1 ) 

TCQ2R=REAL ( CQ2 ) 

TCQ3R=REAL ( CQ3 ) 

TCQ1C=AIMAG(CQ1) 

TCQ2C=AIMAG(CQ2) 

TCQ3C=AIMAG (CQ3 ) 

THOLD ( 1) =DABS (TCQ1R) 

THOLD ( 2 ) =DABS ( TCQ 1C ) 

THOLD ( 3 ) =DABS (TCQ2R) 

THOLD ( 4 ) =DABS (TCQ2C) 

THOLD ( 5 ) =DABS (TCQ3R) 

THOLD ( 6 ) =DABS ( TCQ  3  C ) 

CC - PSI  MOD - REDUCE  COEFFICIENTS  TO  PREVENT  OVERFLOW 

LOGT=0 
DIV=1 . 0 
KOUNT=0 
DO  335  11=1,6 

IF  (THOLD (II) .NE.0.0D0)  THEN 
LOGN=DLOG 1 0 ( THOLD (II)) 

KOUNT=KOUNT+ 1 
LOGT=LOGT+LOGN 


CC 

CC*** 


-•2 _ - 


■  V  V 

■VvV 

•vvSV' 


■\  •*.  •*. 
,*  */  V  V 
A  .  •  :•  *■ 


v.v>>- 

.A.  ir-  n 

•  p.  •  A  A  j 


r  V  v  V 

*w  % 

r  V  V 

v  \*  V 

A  b  •  - 

a;.,. 


•  yAV- 

>’.VA'A 

HTI 


.■«yy 

»*v 


V* «•  < 


A-  3  5 


END  IF 
335  CONTINUE 

IF  (KOUNT.NE. 0)  THEN 
LOGAV=LOGT/KOUNT 
IF  (LOGAV.NE.O)  THEN 
DIV«10.0**LOGAV 
END  IF 
END  IF 

TCQ1R=TCQ1R/DIV 

TCQ1C=TCQ1C/DIV 

TCQ2R=TCQ2R/DIV 

TCQ2C=TCQ2C/DIV 

TCQ3R=TCQ3R/DIV 

TCQ3C=TCQ3C/DIV 

CC - END  OF  PSI  MOD - 

TDSCRR= (TCQ2R*TCQ2R)  -  (TCQ2C*TCQ2C) 

TDSCRR=TDSCRR  -  (4 . 0* ( (TCQ1R*TCQ3R)  -  (TCQ1C*TCQ3C) ) ) 
TDSCRC=2 . 0* (TCQ2R*TCQ2C) 

TDSCRC=TDSCRC  -  (4 . 0* ( (TCQ1R*TCQ3C)  +  (TCQ1C*TCQ3R) ) ) 
CALL  SCROOT  (TDSCRR,TDSCRC, TROOTR, TROOTC) 
TCQ1R=2.0*TCQ1R 
TCQ1C=2.0*TCQ1C 

TDEN=(TCQ1R*TCQ1R)  +  (TCQ1C*TCQ1C) 

TNUMR= ( -TCQ2R)  +  TROOTR 
TNUMC= ( -TCQ2C)  +  TROOTC 

TLSQR=>(  (TNUMR*TCQ1R)  +  (TNUMC*TCQ1C)  )/TDEN 

TLSQC=( (TNUMC*TCQ1R)  -  (TNUMR*TCQ1C) ) /TDEN 

CALL  SCROOT  ( TLSQR , TLSQC , TLR , TLC ) 

A1(I,J)»-434.294481*TLC*2.0 

VI (I , J) = (OMEGA/TLR) /100 . 0 

TNUM= ( -TCQ2R)  -  TROOTR 

TNUMC= (-TCQ2C)  -  TROOTC 

TLSQR=( (TNUMR*TCQ1R)  +  (TNUMC*TCQ1C) ) /TDEN 
TLSQC=( (TNUMC*TCQ1R)  -  (TNUMR*TCQ1C) ) /TDEN 
CALL  SCROOT  (TLSQR, TLSQC, TLR, TLC) 
A2(I,J)=-434.294481*TLC*2 .0 
V2 (I , J) “ (OMEGA/TLR) /100 . 0 
CDEN=CFAC3  -  CFAC2 
C  WRITE  (*,7010)  CDEN 

C  7010  FORMAT  (1X,2E11.4) 

C  WRITE  (*,7011)  DIV 

C  7011  FORMAT  (IX, Ell. 4) 

CTEMPC  -  DSQRT(DIV) 

CTEMPA  =  CDEN/CTEMPC 
CTEMPA  =  CTEMPA/ CTEMPC 
CTEMPB  =  CFAC1/ CTEMPC 
CTEMPB  =  CTEMP  B/ CTEMPC 

CQ1  =  ( (CFAC4*CTEMPA)  -  (CTEMPB*CFAC1) ) / (CSHRB*CTEMPA) 
TCQ 1R=REAL ( CQ 1 ) 

TCQ1C=AIMAG (CQ1) 

CALL  SCROOT  (TCQ1R,TCQ1C,TLR,TLC) 

A3 (I,J)=-434. 29448 1*TLC*2 . 0 
V3 ( I , J) = (OMEGA/TLR) /100 . 0 
CC***  IDENTIFY  TYPE  I,  TYPE  II  DILATATIONAL  WAVES 
IF (VI (I , J) . GE . V2 (I , J) )  GO  TO  350 
EM=V1 ( I , J) 

V1(I,J)-V2(I,J) 

V2 (I , J) =EM 
EM=A1 (I, J) 

A1 ( I , J) =A2 (I,J) 
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cc 

cc 

A2 (I,J)-EM 

v-: 

350 

RETURN 

p 

CC 

CC 

END 

CC 

SUBROUTINE  SCROOT  (TINR,TINC,TOUTR,TOUTC) 

CC 

COMPUTES  DOUBLE  PRECISION  COMPLEX  SQUARE  ROOT 

x 

DOUBLE  PRECISION  TINR,TINC, TRSQ, TCSQ, TSUM,TT,TW, 

CC 

+  TOUTR , TOUTC , DIV 

TO  PREVENT  EXPONENT  OVERFLOW  TAKE  OUT  A  CONSTANT 

DIV* 1.0 

COR* 1.0 

LOGC*0 

*■  ■’ 

< 

IF  (DABS(TINR) . LE . 1 . 0D0 . OR. DABS (TINC) .LE.1.0D0)  i 

Ar 

LOG l=DLOG 1 0 (DABS (TINR) ) 

>  . 

LOG2*DLOG10 (DABS (TINC) ) 

LOGM=LOGl 

IF  (LOG2.LT.LOG1)  LOGM=LOG2 

LOGC=LOGM/2 

DIV*10. ** (LOGC*2) 

COR*10 . **LOGC 

TINR*TINR/DIV 

vl 
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GOTO  1200 


TINC=TINC/DIV 
1200  TRSQ=TINR*TINR 
TCSQ=TINC*TINC 
TSUM=TRSQ+TCSQ 
TW=DSQRT (TRSQ) 

TT=DSQRT( ( (DSQRT(TSUM) ) +TW)/2. 0) 
TW= ( DSQRT ( TCSQ ) )/(2.0*TT) 
SIGNS=-1.0 

IF (TINC.GE .0.0)  SIGNS=1 . 0 
TOUTR=TT 
TOUTC=TW*SIGNS 
IF (TINR.GE. 0.0)  GOTO  1299 
T0UTR=TW 
TOUTCTT*SIGNS 
1299  TINR*TINR*DIV 
TINC*TINC*DIV 
T0UTR*T0UTR*  COR 
TOUTC=TOUTC*COR 
RETURN 
END 


C 

C 

C 

C 

C 

C 

C 

C 

C 


SUBROUTINE  RNDSMP  (ARR,NSIZE,XRND,SMP) 

SUBROUTINE  TO  DETERMINE  VALUES  FOR  THE 
VARIABLES  RHOR,  VOID,  SO,  DECBS ,  DECBE , 
AND  PRATIO  BY  RANDOM  SAMPLING  OF  THEIR 
RESPECTIVE  CUMULATIVE  DISTRIBUTION  CURVES 

DIMENSION  ARR (100 , 2 ) 

CALL  RAND  (XRND) 


RND 


XRND 
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DO  1  I  =  2  ,NSIZE 

IF  (ARR(I,1) .GE.RND)  GOTO  2 

1  CONTINUE 

SMP  -  ARR(NSIZE, 2) 

RETURN 

2  IF  (ARR(I,1) . EQ.RND. AND. I . LT.NSIZE)  THEN 

DO  3  J  =  I+1,NSIZE 

IF  (ARR(J,1) .GT.RND)  GOTO  4 

3  CONTINUE 
J  *  J+l 

4  SMP  =  (ARR(I, 2) +ARR(J-1, 2) )/2 . 0 

ELSE  IF  (ARR(I,1) .EQ.RND. AND. I. EQ.NSIZE)  THEN 
SMP  =  ARR(NSIZE, 2) 

ELSE 

FRACT  =  (RND-ARR(I-1, 1) )/ (ARR(I, 1) -ARR(I-1, 1) ) 
AMOUNT  =  FRACT* (ARR(I , 2) -ARR(I-1, 2) ) 

SMP  =  ARR (1-1,2) + AMOUNT 
END  IF 

RETURN 

END 


SUBROUTINE  RAND  (X) 

A  SUBROUTINE  TO  GENERATE  RANDOM  NUMBERS 
WITH  A  UNIFORM  DISTRIBUTION  OVER  THE  RANGE 
OF  0.0  TO  1.0 


DATA  K,J,M,RM/5701, 3612, 566927, 566927.0/ 

IX  -  INT (X*RM) 

IRAND  -  MOD(J*IX+K,M) 

X  =  (REAL ( IRAND) +0. 5 )/RM 

RETURN 

END 


SUBROTINE  WRIT  ( IFILE , MCRUN , NDEPTH , NFREQ) 

SUBROUTINE  TO  WRITE  VARIOUS  PARAMETERS  OF 
INTEREST  FROM  EACH  MONTE  CARLO  GROUP  OF  RUNS 
AS  WELL  AS  FROM  EACH  INDIVIDUAL  RUN 

COMMON/CONSTS/BLKRR , RHOF , BLKFR , ETA , ALPHA , HOVEMK , 

+  GC1 , GC2 , GC3 , VDH 

COMMON/VDEPTH/EVERY (20) 

COMMON/ VFREQ/ FREQ ( 16 ) 

COMMON/DISTVALS/RHOR, VOID, SO , DECBS , DECBE , PRATIO 
COMMON/ FCTVLS/ BETA , PERM , PS I Z  E 
COMMON/DEPLOOP/EPRESS , G , GPRI , E , EPRI 

COMMON/ FINALS/ A1 (20,16) ,V1(20,16) ,A2(20,16) ,V2(20,16) 
+  A3 (20, 16) ,V3 (20, 16) 


MCPHYS 


IF  (MCRUN.EQ.l)  THEN 

WRITE  (IFILE)  (  EVERY ( I ) ,  I  =  1,20  ) 

WRITE  (IFILE)  (  FREQ (I),  I  =  1,16  ) 

WRITE  (IFILE)  BLKRR,RHOF,BLKFR, ETA, ALPHA 
C  WRITE  (IFILE, 101)  (  EVERY ( I ) ,  I  =  1,20  ) 

C  WRITE  (IFILE, 102)  (  FREQ ( I ) ,  I  =  1,16  ) 

C  WRITE  (IFILE, 103)  BLKRR , RHOF , BLKFR , ETA , ALPHA 

C  101  FORMAT  (20 (F7 . 2 , 2X) ) 

C  102  FORMAT  ( 16 (F7 . 1 , IX) ) 

C  103  FORMAT  (Ell . 5 , IX, F5 . 3 , IX, Ell . 5 , IX, E10 . 4 , IX, F4 . 2) 

END  IF 
C 

WRITE  (IFILE)  MCRUN , RHOR, VOID, SO , DECBS , DECBE , PRATIO 
WRITE  (IFILE)  BETA, PERM, PS IZE , G, GPRI , E, EPRI 
C  WRITE  (IFILE, 200)  MCRUN, RHOR, VOID, SO, DECBS, DECBE, PRATIO 

C  WRITE  (IFILE, 201)  BETA, PERM, PS IZE , G , GPRI , E, EPRI 

DO  10  I  =  1 , NDEPTH 

WRITE  (IFILE)  (  V1(I,J),  J  =  1 , NFREQ  ) 

WRITE  (IFILE)  (  A1(I,J),  J  =  1, NFREQ  ) 

WRITE  (IFILE)  (  V2 (I , J) ,  J  =  1, NFREQ  ) 

WRITE  (IFILE)  (  A2 (I , J) ,  J  =  1, NFREQ  ) 

WRITE  (IFILE)  (  V3 ( I , J) ,  J  =  1 , FREQ  ) 

WRITE  (IFILE)  (  A3 ( I , J) ,  J  =  1, NFREQ  ) 

C  WRITE  (IFILE, 202)  (  V1(I,J),  J  =  1, NFREQ  ) 

C  WRITE  (IFILE, 203)  (  A1(I,J),  J  =  1, NFREQ  ) 

C  WRITE  (IFILE, 202)  (  V2(I,J),  J  =  1, NFREQ  ) 

C  WRITE  (IFILE, 203)  (  A2(I,J),  J  =  1, NFREQ  ) 

C  WRITE  (IFILE, 202)  (  V3(I,J),  J  -  1, NFREQ  ) 

C  WRITE  (IFILE, 203)  (  A3(I,J),  J  =  1, NFREQ  ) 

C  200  FORMAT (15 , IX, 2 (F5 . 3 , IX) , F10 . 3 , IX, 3 (F4 . 3 , IX) ) 

C  201  FORMAT (F6. 4,6(E11.5,1X)) 

C  202  FORMAT (16 (Ell. 5, IX) ) 

C  203  FORMAT (16(E11.5, IX) ) 

10  CONTINUE 
C 
C 

RETURN 
END 

SUBROUTINE  ISREAL(IEND,RNUMB, IERR) 

C  **  ERROR  ERROR 

C  **  CODES  EXPLAINED 

C  ** 

C  ** 

C  ** 

c  ** 

c  ** 

c  ** 

c  ** 

c 

C  **  IMPORTANT:  Subroutine  argument  RNUMB  is  defined  as  REAL*8 . 

C 

character*!  CH(10) 
logical  NUMBER, DECI 
real*4  CHLEFT(IO) ,CHRITE(10) 
real *8  RNUMB 

if  ( IEND  .gt.  10)  then 
IERR  =  6 
RNUMB  =0.0 


1  Contains  at  least  1  non-real  character. 

2  Negative  sign  error. 

3  Imbedded  blank. 

4  Decimal  point  error. 

5  All  blanks  entered. 

6  Argument  "IEND"  is  too  large. 


goto  301 
end  if 

NUMBER  =  .false. 

DECI  =  .false. 

J  -  0 
K  =  0 
SIGN  =  1. 

RNUMB  =  0. 

IERR  =  0 
write  (* , 5) 

format  ('  Enter  Real  number:') 

read  (*,10)  CH(1) ,CH(2) ,CH(3) ,CH(4) ,CH(5) ,CH(6) ,CH(7) ,CH(8) 
CH(9) , CH ( 10) 

format  (lOal) 


do  100 


1, TEND 


if  ( ( CH ( I )  .ne. 

•l*) 

.and. 

( CH ( I )  .ne. 

'2') 

.and. 

& 

( CH ( I )  .ne. 

.3.) 

.and. 

( CH ( I )  .ne. 

.4.) 

.  and. 

& 

( CH ( I )  .ne. 

.51) 

.  and. 

( CH ( I )  .ne. 

•6' ) 

.  and. 

& 

( CH ( I )  .ne. 

.7  .  j 

.  and. 

( CH ( I )  .ne. 

'8') 

.  and. 

& 

( CH ( I )  .ne. 

.9i) 

.and. 

( CH ( I )  .ne. 

'O' ) 

.and. 

& 

& 

( CH ( I )  .ne. 

( CH ( I )  .ne. 
IERR  =  1 
RNUMB  =0.0 
goto  301 

.and. 

then 

( CH ( I )  .ne. 

'  ') 

.  and. 

end  if 


if  ( ( CH ( I ) 

.eq. 

•1') 

.or.  ( CH ( I ) 

.eq. 

'2') 

.or. 

& 

( CH ( I ) 

.eq. 

•3') 

.or.  ( CH ( I ) 

.eq. 

•4') 

.or. 

& 

( CH ( I ) 

.eq. 

'5' ) 

.or.  ( CH ( I ) 

•  eq. 

'  6  ' ) 

.  or. 

& 

( CH ( I ) 

.eq. 

.7. ) 

.or.  ( CH ( I ) 

.eq. 

'8'  ) 

.  or. 

& 

( CH ( I ) 

.eq. 

'9') 

•or.  ( CH ( I ) 

.  eq. 

•0')) 

then 

NUMBER  =  .true, 
if  (.not.  DECI)  then 
J  -  J  +  1 
call  RVALUE (CH (I ) , 
else 

K  =  K  +  1 
call  RVALUE (CH (I) , 
end  if 
goto  100 

else  if  ( CH ( I )  .eq.  '-') 
if  (I  .eq.  IEND)  then 
IERR  =  2 
RNUMB  =0.0 
goto  301 
end  if 

if  (NUMBER)  then 
IERR  =  2 
RNUMB  =0.0 
goto  301 
else 

SIGN  =  -1. 

NUMBER  =  .true, 
goto  100 
end  if 

else  if  ( CH ( I )  .eq.  '  •) 
if  (NUMBER)  then 


CHLEFT(J) ) 


CHRITE(K) ) 


then 


then 


A- 4  0 


£ 

m 

do  80  L  =  I,  IEND 

V\ 

if  ( CH ( 1 )  .ne.  •  •) 
IERR  =  3 

.S 

RNUMB  =0.0 

v 

80 

goto  301 
end  if 
continue 

,v 

else 

i-*  s 

goto  100 
end  if 

s 

else  if  ( CH ( I )  .eq.  '.')  then 

t: 

if  (DECI)  then 

IERR  =  4 

'  jt 

RNUMB  =0.0 

•*  m 

goto  301 

•r , 

else 

DECI  =  .true. 

V 

NUMBER  =  .true. 

goto  100 
end  if 

V, 

end  if 

100 

continue 

V 

t.' 

if  ( (J  .eq.  0)  .and.  (K  .eq.  0) ) 

IERR  =  5 

w  - 

if  (DECI)  IERR  =  4 
if  (SIGN  .eq.  -1.)  IERR  =  2 

RNUMB  =0.0 

goto  301 

end  if 

PLIER  =  1. 

■'/ 

do  200  I  =  J,l,-1 

RNUMB  =  RNUMB  +  (CHLEFT (I)  * 

y. 

PLIER  =  PLIER*  10. 

200 

continue 

m 

VIDER  =  1. 

•/ 

do  300  I  =  1 , K 

VIDER  =  VIDER  *  10. 

RNUMB  =  RNUMB  +  (CHRITE(I)  / 

300 

continue 

RNUMB  =  RNUMB  *  SIGN 

301 

continue 

Ai 

return 

A 

4k 

end 

r* ************************************** 

M 

.v 

SUBROUTINE  RVALUE ( CH , R) 

£ 

character* 1  CH 

real*4  R 

if  (CH  .eq.  '1')  R  =  1. 

MCPHYS 


if  (CH  .eq.  '2')  R  =  2. 
if  (CH  .eq.  '3')  R  =  3. 


A- 41 


PS 
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w 
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Vs.'  V  V 
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'‘VSC^' 
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>-y- 


v'/v 
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■‘-v'v 
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if 

(CH  .eq.  '4') 

R 

= 

4. 

if 

(CH  .eq.  '5') 

R 

ss 

5. 

if 

(CH  .eq.  '6') 

R 

= 

6. 

if 

(CH  .eq.  '7') 

R 

= 

7. 

if 

(CH  .eq.  '8') 

R 

= 

8. 

if 

(CH  .eq.  '9') 

R 

ss 

9. 

if 

(CH  .eq.  '0') 

R 

= 

0. 

return 

end 

c************************************************************ 


SUBROUTINE  GRFDAT (IFILE, INT , MCEND , MCRUN , MIDDEP , MIDFRE ) 
CHARACTER* 1  CHAR (24 ,80) 

COMMON/ FINALS/ A1 (20,16) ,V1(20,16) ,A2(20,16) ,V2(20,16) , 
+  A3(20,16) ,V3(20,16) 

COMMON/VDEPTH/EVERY (20) 

COMMON/ VFREQ/ FREQ ( 16 ) 

DOUBLE  PRECISION  SQSV1 , SQSA1 , SDV1 , SDA1 , SQSV3 , SQSA3 , 

+  SDV3 , SDA3 

IF  (MCRUN  .EQ.  1)  THEN 

I COMP  =  1 
ICOL  =  0 


DO  10  I  =  1,24 

DO  10  J  =  1,80 

CHAR ( I , J )  =  '  .  » 

10  CONTINUE 

SUMV1  =  VI (MIDDEP, MIDFRE) 
SUMA1  =  A1 (MIDDEP, MIDFRE) 
SQSV1  =  VI (MIDDEP, MIDFRE) **2 
SQSA1  =  A1 (MIDDEP, MIDFRE) **2 
RAVGV1  =  SUMV1 
RAVGA1  =  SUMA1 
SDV1  =  0. 

SDA1  =  0. 

A  =  0. 

B  =  0. 

SUMV3  =  V3 (MIDDEP, MIDFRE) 
SUMA3  =  A3 (MIDDEP, MIDFRE) 
SQSV3  =  V3 (MIDDEP, MIDFRE) **2 
SQSA3  =  A3 (MIDDEP, MIDFRE) **2 
RAVGV3  =  SUMV3 
RAVGA3  =  SUMA3 
SDV3  =  0. 

SDA3  *  0. 

C  =  0. 

D  =  0. 

IF  (INT  .EQ.  1)  THEN 
I COMP  =  0 
ICOL  =  1 

IROWA  =  A/.25  +  13 

IROWB  =  B/.25  +  13 

IROWC  =  C/ .25  +  13 

IROWD  =  D/.25  +  13 


p  .v.v.v;^  ’v\" 


3 


<*. 


<. 

A' 

A' 


3 


v 

.v 


MCPHYS 


-v 

\ 


IF  (IROWA.LT. 25. AND. IROWA.GT.O) 
IF  (IROWB.LT. 25. AND. IROWB.GT.O) 
IF  (IROWC.LT. 25. AND. IROWC.GT.O) 
IF  (IROWD.LT. 25. AND. IROWD.GT.O) 
END  IF 


CHAR (IROWA, ICOL) 
CHAR ( IROWB , ICOL) 
CHAR (IROWC, ICOL) 
CHAR (IROWD, ICOL) 


'A' 

'B' 

'C' 

•D' 


ELSE 

ICOMP  =  ICOMP  + 
SUMV1  =  SUMV1  + 
SUMA1  =  SUMA1  + 
SQSV1  =  SQSV1  + 
SQSA1  =  S0SA1  + 


1000 


1001 

50 


1010 


1011 

1012 

90 


VI (MIDDEP,MIDFRE) 
A1(MIDDEP,MIDFRE) 

VI (MIDDEP,MIDFRE) **2 
A1(MIDDEP,MIDFRE) **2 


SQSA1 

RAVGV1  =  SUMV1  /  MCRUN 
RAVGA1  =  SUMA1  /  MCRUN 

SDV1  =  DSQRT (SQSV1  /  MCRUN  -  RAVGV1**2) 
SDA1  =  DSQRT (SQSA1  /  MCRUN  -  RAVGA1**2) 
A  =  DLOGIO (SDV1) 

B  =  DLOGIO (SDA1) 

SUMV3  =  SUMV3  +  V3 (MIDDEP,MIDFRE) 


'123456789 | 1234  56789 | 123456789 | 123456789 |  '  ) 
DO  50  I  =  24, 1,-1 

WRITE  (*,1001)  (CHAR (I , J) , J  =  1 , MCEND) 

FORMAT  (80A1) 

CONTINUE 

IF  (ICOL  .EQ.  MCEND)  THEN 
SCALE  =3.25 
WRITE  (IFILE, 1010) 

FORMAT  ('  123456789 | 123456789 | 123456789 | ' 

'123456789 | 12  3456789 | 12  34  56789 | 123456789 | 123456789 |  '  ) 
DO  90  I  =  24, 1,-1 

SCALE  =  SCALE  -  .25 
WRITE  (IFILE, 1011)  SCALE 
FORMAT  (F5.2) 

WRITE  (IFILE, 1012) (CHAR(I,J) ,J  =  1, MCEND) 
FORMAT  ('  ' , 80A1 ) 

CONTINUE 


SUMA3  =  SUMA3  +  A3 (MIDDEP,MIDFRE) 

» _ — J 

SQSV3  =  SQSV3  +  V3 (MIDDEP,MIDFRE) **2 

SQSA3  =  SQSA3  +  A3 (MIDDEP,MIDFRE) **2 

RAVGV3  =  SUMV3  /  MCRUN 

RAVGA3  =  SUMA3  MCRUN 

SDV3  =  DSQRT (SQSV3  /  MCRUN  -  RAVGV3**2) 

V  V  ' 

SDA3  =  DSQRT (SQSA3  /  MCRUN  -  RAVGA3**2) 

•_  . 

C  =  DLOGIO (SDV3) 

„%V-  .** . 

D  =  DLOGIO (SDA3) 

IF  (INT  .EQ.  ICOMP)  THEN 

ICOMP  =  0 

;>v-> 

ICOL  =  ICOL  +  1 

,‘>v) 

IROWA  =  A/.  25.  +  13 

»  '  ^ 

IROWB  =  B/.25  +  13 

V  V  -> 1 
*'  .  '  ,*• /' 

IROWC  =  C/.25  +  13 

IROWD  =  D/.25  +  13 

IF  (IROWA.LT. 2 5. AND. IROWA.GT.O)  CHAR (IROWA, ICOL) 

=  'A' 

IF  (IROWB.LT. 2 5. AND. IROWB.GT.O)  CHAR (IROWB, ICOL) 

=  'B' 

IF  ( IROWC. LT. 25. AND. IROWC.GT.O)  CHAR ( IROWC , ICOL) 

=  'C' 

's&F 

IF  (IROWD.LT. 2 5. AND. IROWD.GT.O)  CHAR (IROWD, ICOL) 

=  '  D' 

v.  *\  <•.  < 

WRITE  (*,1000) 

■V 

FORMAT  ( ' 123456789 | 123456789 | 123456789 | 123456789 

1  '  , 

/S*  \ 


oooooooooooooo 


SCALE  =  SCALE  -  .25 
WRITE  (IFILE, 1011)  SCALE 

WRITE  (IFILE, 1014)  EVERY (MI DDEP) , FREQ (MI DFRE) 
1014  FORMAT  (//23X, ' DEPTH  =  ' , F7 . 2 , 

+  '  FREQUENCY  =  ' , F9 . 2 ) 

WRITE  (IFILE, 1013)  RAVGV1 , SDV1 , RAVGA1 , SDA1 , 

+  RAVGV3 , S  DV  3 , RAVGA3 , SDA3 


1013 

FORMAT  (/25X,'V1  = 

• ,E9.4,  ' 

+  or  - 

• ,E9 .4/ 

+ 

25X, ' A1  = 

' , E9 . 4 , ' 

+  or  - 

' ,E9. 4/ 

+ 

25X, ' V3  = 

' ,E9.4,  ' 

+  or  - 

• ,E9.4/ 

+ 

25X,'A3  = 

' ,E9.4,  • 

+  or  - 

' ,E9.4) 

END  IF 


END  IF 


END  IF 

WRITE  (*,200)  MCRUN 
WRITE  (*,201)  SUMV1 , SUMA1 
WRITE  (*,201)  SQSV1 , SQSA1 
WRITE  (*,201)  RAVGV1 , RAVGA1 
WRITE  (*,201)  SDV1 , SDA1 
WRITE  (*,202)  A, B 
WRITE  (*,201)  SUMV3 , SUMA3 
WRITE  (*,201)  SQSV3 , SQSA3 
WRITE  (*,201)  RAVGV3 , RAVGA3 
WRITE  (*,201)  SDV3 , SDA3 
WRITE  (*,202)  C, D 

200  FORMAT  (•  MCRUN  -  ',13) 

201  FORMAT  (IX, Ell. 5, IX, Ell. 5) 

202  FORMAT  ( IX, F9 . 6 , IX, F9 . 6) 


RETURN 


,\  «* 
v.v 

V.'j 
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$ LARGE 

PROGRAM  PROFIL 
CHARACTER* 1  YESNO, SCHEME 
CHARACTER* 2  PARNAM 
CHARACTER* 4  PROV,PROF 
CHARACTER*?  TITLE 

CHARACTER* 1 2  SNAME , BNAME1 , BNAME2 , BNAME3 , 

+  FN AME , BN AME  4 , BNAME 5 , BN AME  6 

LOGICAL  FIRST, SOLID, FLUID 

REAL* 8  TMEAN , TSUMSQ , TRTVAR , VALUE , SCORE , SCRLO , SCRHI 
DIMENSION  TMEAN (20, 16,4) , TSUMSQ (20 , 16 , 4 ) , TRTVAR ( 20 , 16 , 4 ) , 

+  WEIGHT (20,16,4) , VALUE ( 20 , 16 , 4 ) , SCORE (2000) , 

+  DEPTH (20) , FREQ ( 16 ) ,V1(20,16) ,A1(20,16) , 

+  V2 (20, 16) ,A2 (20, 16) ,V3 (20, 16) ,A3 (20, 16) , 

+  TNSRMN (20,16,4) , TNSRTP ( 20 , 16 , 4 ) ,TNSRMX(20, 16,4) , 

+  TEMPI (20, 16,4) ,TEMP2 (20, 16,4) , AW(20) , BW( 16) , 

+  CW (4 ) , PARNAM ( 4 ) ,IFPTR(16) ,IDPTR(20) ,THICK(19) 

DIMENSION  WSADEP(20) ,WSAFRQ(16) ,WSAPAR(4) , 

+  WSBDEP(20) ,WSBFRQ(16) ,WSBPAR(4) ,WSGPAR(4) , 

+  WSCDEP (20) , WSCFRQ ( 16 ) ,WSCPAR(4) ,WSFDEP(20) , 

+  WSDDEP (20) , WSDFRQ ( 16 ) ,WSDPAR(4) ,WSEPAR(4) 


+  1  of  file  with  Monte  Carlo  results: 

READ  (*,6001)  PROV 

6001  FORMAT  (A4) 

WRITE  (FNAME , 6002 )  PROV 

6002  FORMAT  ( A4 ,' BIOT . DAT ' ) 

OPEN  ( 2 0,FILE=FNAME,ACCESS=' SEQUENTIAL' , 

+  FORM  =' UNFORMATTED ' ,STATUS=' OLD' ) 

3  WRITE  (*,6005) 


,$) 


6005  FORMAT 
+ 

+ 

+ 

+ 

+ 

+ 


(/ 


Enter  Weighting  Scheme  below'/ 

A  =  All  parameters  weighted  equally  (1.0)'/ 

B  =  Compressional  Speed  &  All  Depths'/ 

weighted  (1.0).  Mid-Frequency  weighted  (1.0)'/ 
C  =  Compressional  Attenuation  &  All  Frequencies'/ 
weighted  (1.0).  Mid-Depth  weighted  (1.0)'/ 

D  =  Compressional  Speed/Attenuation  and  All'/ 

A-  46 


y.\ 
. .  /v 

£ 


C:- 


^ * • 
Mv 
v.v 

•.'.V 
.  -  .  ' 

V ' 


DATA 

WSADEP/6*1. 0,14*0.0/ 

rj  •• 
i\\* 

DATA 

WSAFRQ/7*1. 0, 9*0.0/ 

k’/' 

DATA 

WS APAR/ 4*1.0/ 

v/L 

%  -*• 

>  **-» 

DATA 

WSBDEP/6*1. 0, 14*0. 0/ 

a# 

DATA 

WSBFRQ/3*0. 0,1. 0,12*0.0/ 

.\\\ 

DATA 

WSBPAR/1.0, 3*0.0/ 

■/* 

'  m+ 

* 

DATA 

WSCDEP/3*0. 0,1. 0,16*0.0/ 

•* 

*»  /* 

r-s 

DATA 

WSCFRQ/7*1.0, 9*0.0/ 

DATA 

WSCPAR/0. 0,1.0, 2*0.0/ 

'  ^ 

DATA 

WSDDEP/6*1.0, 14*0.0/ 

n 

;  * 

[V 

DATA 

WSDFRQ/7*1.0, 9*0.0/ 

«  K 

Cv-' 

DATA 

WS  DPAR/ 2*1.0, 2*0.0/ 

\y.- 

DATA 

WSEPAR/0 . 0,0. 0,1. 0,0.0/ 

•\V 

DATA 

WSFDEP/2*0. 0,1. 0,3 *0.0, 14*0.0/ 

DATA 

WSGPAR/3*0 .0,1.0/ 

*  v 

DATA 

PARNAM/ ' VI ' , ' A1 ' , ' V3 ' , ' A3 ' / 

*. 

ft 

FIRST  =  .TRUE. 

* 

WRITE  (*,6000) 

$ 

*81111111 

6000  FORMAT  (/ '  Enter  4  character  province  code'/ 

*  “* 

•  S 

:S  ^ 

-  i:  «! 

» * 

"  % 
,  £:■ 


'  V  V 


PROFIL 

Depths  and  Frequencies  weighted  (1.0)'/ 

E  =  Shear  Speed  and  All  Depths  weighted  (1.0).'/ 
Mid-Frequency  weighted  (1.0).'/ 

F  =  Compress ional  Attenuation  &  All  Frequencies'/ 
weighted  (1.0).  5m.  Depth  weighted  (1.0).'/ 

G  =  Shear  Attenuation  and  All  Frequencies'/ 

weighted  (1.0).  5m.  Depth  weighted  (1.0).'/ 
Enter  Scheme:  ',$) 


READ  (*,6105)  SCHEME 
6105  FORMAT  (Al) 


IF  (SCHEME  .EQ. 
SCHEME  =  'A' 

•A' 

•  OR. 

SCHEME  .EQ.  'a')  THEN 

n 

.V. 

ELSE  IF  (SCHEME 
SCHEME  =  'B' 

•  EQ. 

•B' 

•  OR.  SCHEME  .EQ.  ' b  * ) 

THEN 

ELSE  IF  (SCHEME 
SCHEME  =  'C' 

.EQ. 

•c* 

.OR.  SCHEME  .EQ.  'c') 

THEN 

,V 

ELSE  IF  (SCHEME 
SCHEME  =  'D' 

.EQ. 

'D' 

.OR.  SCHEME  .EQ.  'd') 

THEN 

IS 

ELSE  IF  (SCHEME 
SCHEME  =  'E' 

•  EQ. 

'E' 

.OR.  SCHEME  .EQ.  'e') 

THEN 

tl 

ELSE  IF  (SCHEME 
SCHEME  =  'F' 

•  EQ. 

•F' 

.OR.  SCHEME  .EQ.  ' f ' ) 

THEN 

A 

ELSE  IF  (SCHEME 

•  EQ. 

'G' 

•OR.  SCHEME  .EQ.  'g') 

THEN 

SCHEME  =  'G' 
ELSE 

GOTO  3 
END  IF 


WRITE  (SNAME, 6006)  PRO V, SCHEME 

6006  FORMAT  ( A4 , ' WS ' , Al , ' . PRO ' ) 

OPEN  ( 10 , FILE=SNAME , ACCESS* 'SEQUENTIAL ' , 
+  FORM* ' FORMATTED ' , STATUS* ' NEW ' ) 

WRITE  (TITLE, 6007)  PROV, SCHEME 

6007  FORMAT  (A4,'WS',A1) 


5  WRITE  (*,6010) 

6010  FORMAT  (/ '  Enter  the  number  of  Monte  Carlo  steps'/ 
+  '  performed  to  produce  the  above'/ 

+  '  mentioned  results:  ',$) 

READ  (*, 6015, err=5)  MCLOOP 
6015  FORMAT  (14) 


DO  10  I  =  1,20 
DO  10  J  *  1,16 
DO  10  K  =  1,4 

TMEAN (I , J, K)  =0.0 
TSUMSQ ( I , J , K)  =  0.0 
TRTVAR ( I , J , K)  =  0.0 
VALUE ( I , J, K)  =  0.0 

10  CONTINUE 

READ  (20)  (DEPTH (I) ,  I  =  1,20) 

READ  (20)  (FREQ (I) ,  I  =  1,16) 

READ  (20)  BLKRR,RHOF, BLKFR, ETA, ALPHA 

DO  19  I  =  2,20 

IF  ( DEPTH ( I )  .LE.  0.0)  GOTO  20 
19  CONTINUE 
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20  NDEPTH  =1-1 


DO  29  I  =  2,16 

IF  (FREQ (I)  .LE.  0.0)  GOTO  30 

29  CONTINUE 
1  =  1  +  1 

30  NFREQ  =1-1 

SOLID  =  .FALSE. 

FLUID  =  .FALSE. 


1000  WRITE  (*,6110) 

611 J  FORMAT  (/ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

READ  (*,6105)  YESNO 
IF  (YESNO  .EQ.  'A' 

SOLID  =  .TRUE. 

ELSE  IF  (YESNO  .EQ. 

FLUID  =  .TRUE. 

ELSE  IF  (YESNO  .EQ.  'C' 
SOLID  =  .TRUE. 

FLUID  =  .TRUE. 

ELSE  IF  (YESNO  .EQ.  'D' 
GOTO  31 
ELSE 

GOTO  1000 
END  IF 


',$) 


YESNO  .EQ.  'a') 

THEN 

.OR. 

YESNO  .EQ. 

'b' ) 

THEN 

.OR. 

YESNO  .EQ. 

•c' ) 

THEN 

.OR. 

YESNO  .EQ. 

'd' ) 

THEN 

Choose  from  below: */ 

A  =  Create  Multilayer  (Solid)  Model'/ 
input  files  (Min,  Max,  Typ) '/ 

B  =  Create  Reflec  (Fluid)  Model'/ 
input  files  (Min,  Max,  Typ) '/ 

C  =  Create  Both  (Solid  &  Fluid)  Model'/ 
input  files  (Min,  Max,  Typ) '/ 

D  =  None  of  the  above. '/ 

Enter  choice: 

.OR. 

'B' 


WRITE  (*,6115) 

6115  FORMAT  (/ '  Choose  Frequencies  for  Bottom  Loss'/ 

+  '  profiles  from  the  following  list:') 

LBFREQ  =  0 

DO  1010  I  =  1, NFREQ 

WRITE  (*,6120)  FREQ ( I ) 

6120  FORMAT  (1X,G9.3,'  Hz.  [Y/N]?  ',$) 

READ  (*,6105)  YESNO 

IF  (YESNO  .EQ.  'N'  .OR.  YESNO  .EQ.  'n')  GOTO  1010 
LBFREQ  =  LBFREQ  +  1 
I FPTR( LBFREQ)  =  I 
1010  CONTINUE 

1014  WRITE  (*,6121) 

6121  FORMAT  (/ '  Choose  Depths  (All  or  One)  for  Bottom'/ 

+  '  Loss  profiles  from  the  following  list: '/ 

+  '  0  =  All  Depths') 

DO  1015  I  =  1, NDEPTH 

WRITE  (*,6122)  I,  DEPTH (I) 

6122  FORMAT  (IX, 12,'  =  ',G9.3) 

IDPTR(I)  =  I 

1015  CONTINUE 
WRITE  (*,6123) 

6123  FORMAT  ('  Enter  choice:  ',$) 


v. 


»TT»T?  V  V’AV.’ 


PROFIL 

READ  (*,6124)  IVAL 
6124  FORMAT  (12) 

IF  (IVAL  .LT.  0  .OR.  IVAL  .GT.  NDEPTH)  GOTO  1014 
IF  (IVAL  .NE.  0)  THEN 
NUMLAY  =  1 

THICK (NUMLAY)  =  DEPTH (NDEPTH) 

IDPTR (NUMLAY)  =  IVAL 
ELSE 

NUMLAY  =  NDEPTH  -  1 
PRETHK  =  0.0 
DO  1016  I  -  1, NUMLAY 

TEMP  =  (DEPTH(I+1)  -  DEPTH ( I ) )  /  2.0 
THICK (I)  =  TEMP  +  PRETHK 
PRETHK  =  TEMP 
1016  CONTINUE 
END  IF 

ANGMIN  =  1.0 


.'.v.v , 

m 
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k  *  • m  * 

*  ;.*  v 


vy.v 
> v* 


ANGMAX  =90.0 
ANGINC  =1.0 

IF  (SOLID)  THEN 

WRITE  ( BNAME 1,6160)  PRO V, SCHEME 
6160  FORMAT  ( 'S' ,A4, 'WS' ,A1, ' .MIN' ) 

OPEN  (30 , FILE=BNAME1 , ACCESS= ' SEQUENTIAL' , 
+  FORM= ' FORMATTED ' , STATUS = ' NEW ' ) 

WRITE  (BNAME2 , 6165)  PROV, SCHEME 
6165  FORMAT  ( ' S • , A4 , ' WS • , Al , ' .TYP' ) 

OPEN  ( 40, FILE=BNAME2,ACCESS=' SEQUENTIAL' , 
+  FORM= ' FORMATTED ' , STATUS* ' NEW ' ) 

WRITE  (BNAME3 ,6170)  PROV, SCHEME 
6170  FORMAT  ( ' S ' , A4 , ' WS ' , Al , ' .MAX ' ) 

OPEN  ( 5 0 , FILE=BNAME3 , ACCESS* ' SEQUENTIAL' , 
+  FORM* ' FORMATTED ' , STATUS* ' NEW ' ) 

END  IF 

IF  (FLUID)  THEN 

WRITE  (BNAME4 ,6175)  PROV, SCHEME 
6175  FORMAT  ( ' F ' , A4 , ' WS ' , Al , ' .MIN' ) 

OPEN  ( 60 , FILE=BNAME4 , ACCESS* ' SEQUENTIAL ' , 
+  FORM* ' FORMATTED ', STATUS* 'NEW' ) 

WRITE  (BNAME5 ,6180)  PROV, SCHEME 
6180  FORMAT  ( ' F ' , A4 , ' WS ' , Al, ' .TYP ' ) 

OPEN  ( 7 0,FILE=BNAME5, ACCESS* 'SEQUENTIAL* , 
+  FORM*' FORMATTED' , STATUS* ' NEW' ) 
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WRITE  (BNAME 6, 6185)  PROV, SCHEME 
6185  FORMAT  ( ' F ' , A4 , ' WS ' , Al , ' .MAX ' ) 

OPEN  (80, FILE=BNAME6, ACCESS* 'SEQUENTIAL' , 
+  FORM* ' FORMATTED ' , STATUS* ' NEW ' ) 

END  IF 


31  DO  32 

I  *  1, NDEPTH 

— 

IF 

(SCHEME  .EQ. 

•A’) 

AW ( I)  *  WSADEP ( I ) 

IF 

(SCHEME  .EQ. 

'  B  *  ) 

AW(I)  *  WSBDEP ( I ) 

IF 

(SCHEME  .EQ. 

'  C ' ) 

AW ( I )  =  WSCDEP ( I ) 

IF 

(SCHEME  .EQ. 

'D' ) 

AW ( I )  =  WSDDEP(I) 

IF 

(SCHEME  .EQ. 

'  E ' ) 

AW(I)  =  WSADEP ( I) 

A- 4  9 
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IF  (SCHEME  .EQ.  'F')  AW(I)  =  WSFDEP(I) 

IF  (SCHEMEM.EQ.  'G')  AW(I)  =  WSFDEP(I) 

32  CONTINUE 

DO  34  I  =  1,NFREQ 

IF  (SCHEME  .EQ.  'A')  BW(I)  =  WSAFRQ(I) 

IF  (SCHEME  .EQ.  'B')  BW(I)  =  WSBFRQ(I) 

IF  (SCHEME  .EQ.  'C')  BW(I)  =  WSCFRQ(I) 

IF  (SCHEME  .EQ.  'D')  BW(I)  =  WSDFRQ(I) 

IF  (SCHEME  .EQ.  * E ' )  BW(I)  =  WSBFRQ(I) 

IF  (SCHEME  .EQ.  'F')  BW(I)  =  WSAFRQ(I) 

IF  (SCHEME  .EQ.  'G*)  BW(I)  =  WSAFRQ(I) 

34  CONTINUE 

DO  36  I  =  1,4 

IF  (SCHEME  .EQ.  'A')  CW(I)  =  WSAPAR(I) 

IF  (SCHEME  .EQ.  *  B ' )  CW(I)  =  WSBPAR(I) 

IF  (SCHEME  .EQ.  'C')  CW(I)  =  WSCPAR(I) 

IF  (SCHEME  .EQ.  'D')  CW(I)  =  WSDPAR(I) 

IF  (SCHEME  .EQ.  'E*)  CW(I)  =  WSEPAR(I) 

IF  (SCHEME  .EQ.  *F')  CW(I)  =  WSCPAR(I) 

IF  (SCHEME  .EQ.  'G')  CW(I)  =  WSGPAR(I) 

36  CONTINUE 

DO  38  I  =  1 , NDEPTH 
DO  38  J  =  1 , NFREQ 
DO  38  K  =  1,4 

WEIGHT ( I , J , K)  =  AW(I)  *  BW(J)  *  CW(K) 

38  CONTINUE 

DO  100  M  *  1 , MCLOOP 
SCORE (M)  =0.0 

READ  (20)  MCRUN , RHOR, VOID, SO, DECBS , DECBE , PRATIO 
READ  (20)  BETA, PERM, PSIZE,G,GPRI,E,EPRI 
DO  40  I  =  1, NDEPTH 

READ  (20)  (V1(I,J),  J  =  1, NFREQ) 

READ  (20)  (A1(I, J) ,  J  =  1, NFREQ) 

READ  (20)  (V2 (I , J) ,  J  =  1, NFREQ) 

READ  (20)  (A2 (I , J) ,  J  =  1, NFREQ) 

READ  (20)  (V3(I,J),  J  =  1, NFREQ) 

READ  (20)  (A3 (I , J) ,  J  =  1, NFREQ) 

DO  40  K  =  1, NFREQ 

TMEAN ( I , K , 1 )  =  TMEAN (I , K, 1)  +  V1(I,K) 

TME AN ( I , K , 2 )  =  TMEAN ( I , K , 2 )  +  A1(I,K) 

TMEAN ( I , K , 3 )  =  TMEAN ( I , K , 3 )  +  V3(I,K) 

TMEAN (I, K, 4)  =  TMEAN ( I , K , 4 )  +  A3(I,K) 
TSUMSQ ( I , K , 1 )  =  TSUMSQ (I , K, 1)  +  V1(I,K)**2 
TSUMSQ ( I , K, 2 )  =  TSUMSQ (I , K, 2 )  +  A1(I,K)**2 
TSUMSQ (I , K, 3 )  =  TSUMSQ ( I , K, 3 )  +  V3(I,K)**2 
TSUMSQ ( I , K, 4 )  =  TSUMSQ (I, K, 4)  +A3(I,K)**2 
40  CONTINUE 

100  CONTINUE 

DO  110  I  =  1, NDEPTH 
DO  110  J  =  1, NFREQ 
DO  110  K  =  1,4 

TMEAN (I, J,K)  =  TMEAN ( I , J, K)  /  MCLOOP 
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VALUE ( I ,  J ,  K)  =  TMEAN(I , J,K) 

TRTVAR ( I ,  J ,  K)  =  DSQRT ( (TSUMSQ ( I , J , K)  -  MCLOOP  * 

+  (TMEAN (I, J,K)  **  2))/  (MCLOOP  -  1) ) 

110  CONTINUE 

6023  FORMAT  ('  ') 

115  CONTINUE 

SCRLO  =  1.0E35 

SCRHI  =  -1.0E35 

REWIND  (20) 

READ  (20)  (DEPTH (I) ,  I  =  1,20) 

READ  (20)  (FREQ (I) ,  I  =  1,16) 

READ  (20)  BLKRR,RHOF,BLKFR, ETA, ALPHA 

DO  200  M  =  1, MCLOOP 

READ  (20)  MCRUN, RH0R, VOID, SO , DECBS , DECBE, PRATIO 
READ  (20)  BETA, PERM, PSIZE,G,GPRI,E,EPRI 
DO  140  I  =  1 , NDEPTH 

READ  (20)  (VI (I , K) ,  K  =  1,NFREQ) 

READ  (20)  (A1 (I , K) ,  K  =  1,NFREQ) 

READ  (20)  (V2(I,K),  K  =  1, NFREQ) 

READ  (20)  ( A2 ( I , K) ,  K  =  1, NFREQ) 

READ  (20)  (V3 ( I , K) ,  K  =  1, NFREQ) 

READ  (20)  (A3 (I,K) ,  K  =  1, NFREQ) 

DO  140  J  -  1, NFREQ 

IF  (TRTVAR(I, J, 1)  .NE.  0.0) 

+  SCORE (M)  =  (DABS (VI (I , J)  -  VALUE ( I , J , 1 ) )  / 

+  TRTVAR ( I , J , 1 ) )  *  WEIGHT ( I , J , 1 )  + 

+  SCORE (M) 

IF  ( TRTVAR ( I , J , 2 )  .NE.  0.0) 

+  SCORE (M)  =  (DABS (A1(I, J)  -  VALUE (I , J , 2) )  / 

+  TRTVAR ( I , J , 2 ) )  *  WEIGHT(I,J,2)  + 

+  SCORE (M) 

IF  (TRTVAR ( I , J , 3 )  .NE.  0.0) 

+  SCORE (M)  =  (DABS (V3 (I , J)  -  VALUE ( I , J , 3 ) )  / 

+  TRTVAR ( I , J , 3 ) )  *  WEIGHT ( I, J, 3)  + 

+  SCORE (M) 

IF  (TRTVAR ( I , J , 4 )  .NE.  0.0) 

+  SCORE (M)  =  (DABS (A3 (I, J)  -  VALUE ( I , J , 4 ) )  / 

+  TRTVAR ( I , J , 4 ) )  *  WEIGHT(I,J,4)  + 

+  SCORE (M) 

140  CONTINUE 

IF  (SCORE (M)  .LT.  SCRLO)  THEN 
SCRLO  =  SCORE (M) 

MVALLO  *  M 
DO  150  I  =  1, NDEPTH 
DO  150  J  =  1, NFREQ 

RHOR1  =  RHOR 

BETA1  =  BETA 

TEMPI (I , J, 1)  =  V1(I,J) 

TEMPI (I , J, 2)  =  A1 ( I , J) 

TEMPI ( I , J , 3 )  =  V3 (I , J) 

TEMPI (I ,J, 4 )  =  A3 (I , J) 

150  CONTINUE 

END  IF 
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IF  (SCORE (M)  .GT.  SCRHI)  THEN 
SCRHI  =  SCORE (M) 

MVALHI  =  M 
DO  160  I  =  1,NDEPTH 
DO  160  J  =  1 , NFREQ 
RHOR2  =  RHOR 
BETA2  =  BETA 
TEMP2 (I , J, 1)  =  V1(I,J) 
TEMP2 ( I ,  J ,  2 )  =  A1 ( I ,  J) 
TEMP2 ( I ,  J ,  3 )  =  V3 ( I , J) 
TEMP2 ( I ,  J ,  4 )  =  A3 (I , J) 
CONTINUE 
END  IF 


SCORE (M)  =0.0 

200  CONTINUE 

DO  240  I  =  1 , NDEPTH 
DO  240  J  =  1, NFREQ 
DO  240  K  =  1,4 

IF  (FIRST)  THEN 

TNSRTP ( I , J , K)  =  TEMPI (I , J , K) 

TNSRMX ( I , J , K)  =  TEMP2 ( I , J , K) 

VALUE ( I , J , K)  =  DBLE (TNSRMX (I, J,K) ) 


RHORTP  =  RHOR1 
BETATP  =  BETA1 
RHORMX  =  RHOR2 
BETAMX  =  BETA2 
SCRTYP  =  SNGL(SCRLO) 

MTPVAL  =  MVALLO 
SCRMAX  =  SNGL( SCRHI) 

MAXVAL  =  MVALHI 
ELSE 

TNSRMN ( I , J , K)  =  TEMP2 ( I , J , K) 
RHORMN  =  RHOR2 
BETAMN  =  BETA2 
SCRMIN  =  SNGL( SCRHI) 

MINVAL  =  MVALHI 
MSAME  =  MVALLO 
SCRNEW  =  SNGL(SCRLO) 

END  IF 


240  CONTINUE 

IF  (FIRST)  THEN 
FIRST  =  .FALSE. 

GOTO  115 
END  IF 

WSPEED  =  SQRT ( BLKFR/RHOF)  /  100.0 

DENSMN  =  RHOF  *  BETAMN  +  RHORMN  *  (1  -  BETAMN) 

DENSTP  =  RHOF  *  BETATP  +  RHORTP  *  (1  -  BETATP) 

DENSMX  =  RHOF  *  BETAMX  +  RHORMX  *  (1  -  BETAMX) 

IF  (SOLID)  THEN 

WRITE  (30,6500)  ANGMIN , ANGMAX , ANGINC 
WRITE  (40,6500)  ANGMIN, ANGMAX, ANGINC 
WRITE  (50,6500)  ANGMIN, ANGMAX, ANGINC 


6500 


FORMAT  (3F10.1) 
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PROFIL 


WRITE  (30,6505)  WSPEED, RHOF , LBFREQ 
WRITE  (40,6505)  WSPEED, RHOF , LBFREQ 
WRITE  (50,6505)  WSPEED, RHOF, LBFREQ 
6505  FORMAT  (F10. 1, F10. 3 , 14) 

DO  245  I  -  1, LBFREQ 

TF  -  FREQ(IFPTR(I) )/1000. 0 
WRITE  (30,6510)  TF , NUMLAY 
WRITE  (40,6510)  TF, NUMLAY 
WRITE  (50,6510)  TF, NUMLAY 
6510  FORMAT  (F10.4,I4) 

WRITE  (30,6515)  TNSRMN (NDEPTH, IFPTR(I) , 1) , 

TNSRMN (NDEPTH , IFPTR (I) ,3) , DENSMN 
WRITE  (40,6515)  TNSRTP (NDEPTH, IFPTR(I) , 1) , 

TNSRTP (NDEPTH , IFPTR (I) ,3) , DENSTP 
WRITE  (50,6515)  TNSRMX (NDEPTH, IFPTR (I) , 1) , 

TNSRMX( NDEPTH, IFPTR (I) ,3) , DENSMX 
FORMAT  (2F10 . 1, F10 . 3 ) 

CAVMN  =  TNSRMN (NDEPTH, IFPTR (I) ,2)  /  (FREQ (IFPTR(I) ) /1000 . ) 
SAVMN  =  TNSRMN (NDEPTH, IFPTR (I) ,4)  /  (FREQ (IFPTR (I) )/1000. ) 
CAVTP  =  TNSRTP (NDEPTH, IFPTR (I) ,2)  /  (FREQ (IFPTR (I) )/1000. ) 
SAVTP  =  TNSRTP (NDEPTH, IFPTR(I) ,4)  /  (FREQ (IFPTR (I) )/1000 . ) 
CAVMX  =  TNSRMX( NDEPTH, IFPTR (I) ,2)  /  (FREQ (IFPTR (I) ) /1000 . ) 
SAVMX  =  TNSRMX (NDEPTH, IFPTR (I) ,4)  /  (FREQ ( IFPTR ( I) ) /1000 . ) 
WRITE  (30,6520)  CAVMN, SAVMN 
WRITE  (40,6520)  CAVTP, SAVTP 
WRITE  (50,6520)  CAVMX, SAVMX 
6520  FORMAT  (2F10.4) 

DO  245  J  =  1, NUMLAY 

WRITE  (30,6525)  TNSRMN (IDPTR(J) , IFPTR (I) , 1) , 

TNSRMN (IDPTR(J) , IFPTR (I) ,3) , DENSMN 
WRITE  (40,6525)  TNSRTP ( IDPTR (J) , IFPTR ( I ), 1) , 

TNSRTP (IDPTR(J) , IFPTR (I) ,3) , DENSTP 
WRITE  (50,6525)  TNSRMX ( IDPTR (J) , IFPTR(I) , 1) , 

TNSRMX ( IDPTR (J) ,IFPTR(I) ,3) , DENSMX 
FORMAT  (F10.1,F10.2,F10.3) 

CAVMN=TNSRMN( IDPTR (J) ,IFPTR(I) , 2 ) / ( FREQ (IFPTR ( I ) )/1000.) 
SAVMN=TNSRMN( IDPTR (J) , IFPTR (I) ,4)/ (FREQ(IFPTR(I) ) /1000. ) 
CAVTP=TNSRTP( IDPTR (J) ,IFPTR(I) , 2)/ (FREQ (IFPTR (I) )/1000. ) 
SAVTP=TNSRTP( IDPTR (J) ,IFPTR(I) , 4)/ (FREQ (IFPTR (I) )/1000. ) 
CAVMX=TNSRMX( IDPTR (J) , IFPTR (I) , 2 ) / (FREQ ( IFPTR ( I) ) /1000 . ) 
SAVMX=TNSRMX( IDPTR (J) ,IFPTR(I) , 4 ) / (FREQ ( IFPTR ( I) )/1000. ) 
TTH  -  THICK (J)  *  100.0 
WRITE  (30,6530)  CAVMN, SAVMN, TTH 
WRITE  (40,6530)  CAVTP , SAVTP , TTH 
WRITE  (50,6530)  CAVMX , SAVMX , TTH 
6530  FORMAT  ( 2F10 . 4 , F10 . 1) 

245  CONTINUE 

CLOSE  (30) 

CLOSE  (40) 

CLOSE  (50) 

END  IF 


+ 

+ 

+ 
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IF  (FLUID)  THEN 


WRITE  (60,6124)  LBFREQ 
WRITE  (70,6124)  LBFREQ 
WRITE  (80,6124)  LBFREQ 
IFILE  =  1 

WRITE  (60,6600)  IFILE, BNAME4 
WRITE  (70,6600)  IFILE, BNAME5 
WRITE  (80,6600)  IFILE, BNAME6 
6600  FORMAT  (I5,A50) 

DO  2900  I  =  1, LBFREQ 
KEY1  =  1 
KEY 2  =  0 
KEY 3  =  0 
KEY PAR  =  3 

WRITE  (60,6605)  KEY1 , KEY2 , KEY 3 , KEYPAR 
WRITE  (70,6605)  KEY1 , KEY2 , KEY3 , KEYPAR 
WRITE  (80,6605)  KEY 1, KEY 2 , KEY 3 , KEYPAR 
6605  FORMAT  (415) 

NCASE  =  1 
NANG  =  90 
BANG  *  1.0 
DINC  =1.0 
FLO  =  FREQ (I) 

FHI  =  FREQ ( I ) 

NPTF  =  1 

WRITE  (60,6610)  NCASE, NANG, BANG, DINC, FLO, FHI , NPTF 
WRITE  (70,6610)  NCASE, NANG, BANG, DINC , FLO , FHI , NPTF 
WRITE  (80,6610)  NCASE, NANG, BANG, DINC, FLO , FHI , NPTF 
6610  FORMAT  ( 215 , 2F10 . 2 , 2F10 . 2 , 15) 

HC3  =  -1.0 
HALF1  =0.0 
HALF 3  =0.0 

CAMN  =  TNSRMN (NDEPTH, IFPTR(I) , 2 ) / (FREQ (IFPTR (I) ) /1000 . ) 
CATP  =  TNSRTP (NDEPTH, IFPTR(I) ,2)/ (FREQ (IFPTR (I) ) /1000. ) 
CAMX  =  TNSRMX (NDEPTH, IFPTR (I) , 2 ) / (FREQ (IFPTR (I) ) /1000 . ) 
WRITE  (60,6615)  WSPEED, TNSRMN (NDEPTH, IFPTR(I) , 1) , 

+  HC3 , RHOF, DENSMN , HALF1 , CAMN, HALF3 

WRITE  (70,6615)  WSPEED, TNSRTP (NDEPTH, IFPTR( I) , 1) , 

+  HC3, RHOF, DENSTP,HALF1, CATP, HALF3 

WRITE  (80,6615)  WSPEED, TNSRMX (NDEPTH, IFPTR(I) ,1) , 

+  HC3, RHOF, DENSMX,HALF1, CAMX, HALF3 

6615  FORMAT  (3F10 . 2 , 2F10 . 3 , 2F10 . 8 , F10 . 7) 

WRITE  (60,6620)  NUMLAY 

WRITE  (70,6620)  NUMLAY 

WRITE  (80,6620)  NUMLAY 

6620  FORMAT  (15) 

DO  2899  J  =  1, NUMLAY 

WRITE  (60,6625)  TNSRMN ( IDPTR (J) , IFPTR (I) , 1) , DENSMN 
WRITE  (70,6625)  TNSRTP ( IDPTR (J) , IFPTR (I) , 1) , DENSTP 
WRITE  (80,6625)  TNSRMX ( IDPTR (J) , IFPTR(I) , 1) , DENSMX 
CAMN  =  TNSRMN ( IDPTR (J) , IFPTR (I) ,2)  / 

+  (FREQ (IFPTR ( I ) ) /1000 . ) 

CATP  =  TNSRTP ( IDPTR (J) ,IFPTR(I) ,2)  / 

+  (FREQ(IFPTR(I) )/1000. ) 

CAMX  =  TNSRMX ( IDPTR (J) ,IFPTR(I) ,2)  / 

+  (FREQ (IFPTR (I) )/1000. ) 


PROFIL 


6625 

2899 

WRITE 

WRITE 

WRITE 

FORMAT 

CONTINUE 

(60.6625) 

(70.6625) 

(80.6625) 
(2F10.4) 

CAMN, THICK (J) 
CATP, THICK (J) 
CAMX, THICK (J) 

2900 

CONTINUE 

END  IF 

WRITE  (10,6029)  TITLE 

6029  FORMAT  (A7,'  PROFILE') 

WRITE  (10,6023) 

WRITE  (10,6030)  MTPVAL,SCRTYP,MAXVAL, SCRMAX 

6030  FORMAT  ('SCORE  ',14,'  TYPICAL  =  ',G12.4, 

+  '  SCORE  ',14,'  MAXIMUM  =  ',G12.4) 

WRITE  (10,6031)  MSAME , SCRNEW , MINVAL , SCRMIN 

6031  FORMAT  ('SCORE  ',14,'  NEW  =  • ,G12.4, 

+  '  SCORE  ',14,'  MINIMUM  =  ',G12.4) 

WRITE  (10,6095)  NDEPTH 
6095  FORMAT  (12) 

DO  250  I  =  1, NDEPTH 

WRITE  (10,6075)  DEPTH ( I ) ,  AW(I) 

6075  FORMAT  ('DEPTH  ',F9.3,'  WEIGHT  FACTOR  =  ',F4.2) 

250  CONTINUE 

WRITE  (10,6095)  NFREQ 
DO  260  I  =  1, NFREQ 

WRITE  (10,6080)  FREQ (I) ,  BW(I) 

6080  FORMAT  ('FREQUENCY  ' , F9 . 3 , •  WEIGHT  FACTOR  =  ',F4.2) 

260  CONTINUE 

WRITE  (10,6023) 

DO  270  1=1,4 

WRITE  (10,6085)  PARNAM ( I ) ,  CW(I) 

6085  FORMAT  ('PARAMETER  ',A2,'  WEIGHT  FACTOR  =  ',F4.2) 

270  CONTINUE 

WRITE  (10,6023) 

WRITE  (10,6090) 

6090  FORMAT  (6X, ' DEPTH(I) ', 4X, ' FREQUENCY (J) ', 4X, ' PARAMETER(K) ' , 
+  7X, 'TMEAN' ,11X, ' TRTVAR ' ,10X, 'TNSRMN ' ,10X, ' TNSRTP '  , 

+  10X, 'TNSRMX' ) 

WRITE  (10,6023) 

DO  300  I  =  1, NDEPTH 
DO  300  J  =  1, NFREQ 
DO  300  K  =  1,4 

WRITE  (10,6035)  I , J , K, TMEAN (I , J , K) , TRTVAR ( I , J , K) 
+  TNSRMN ( I ,  J ,  K) , TNSRTP ( I , J , K) , 

+  TNSRMX ( I , J , K) 

6035  FORMAT  (9X, 12 , 12X, 12 , 14X, II , 10X, 5 (E13 . 7 , 3X) ) 

300  CONTINUE 

CLOSE  (20) 

CLOSE  (10) 
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First  Letter 
(seafloor  depth) 


Letter 

A 

B 

C 

Z 


Depth 

Interval  (m) 
0-200 

201-2000 

2001-4000 

All  depths 


Second  Letter 
(grain  size  type) 


Letter 

Name(s) 

Class 
Number  (s) 

A 

Sand 

1 

B 

Silty  sand 

2 

Muddy  sand 

3 

Clayey  sand 

4 

C 

Sandy  silt 

5 

Sandy  mud 

6 

Sandy  clay 

7 

D 

Silt 

8 

E 

Mud 

9 

Clay 

10 

Z 


All  size 
types 


1  thru  10 


PROVINCE:  AAZZ 


APPENDIX  B 


r.'  *v  *  v  »  j« 


'J  -  - 

T  T  V 

<  ^  -1- 

1  % 


./.V. 


.V  . 

■'■Vi 
a\.vV- . 

r  . 


%$ 

.  *  *'* 

»  ^ 

If 

a’a  a- 

c-:c-: 


y 


•j 

-j 


:-vv*i 

k _ 


\W 


£ 
sr:f!B 


PROVINCE:  ABZZ 

KKtOULNCY:  1600.0  Hi.  GKA.iING  ANGLE:  S.O 


PROVINCE:  ACZZ 


PROVINCE:  ADZZ 

f-RLOMLNCY:  100.(5  Hz.  GRAZING  ANGLE: 


PROVINCE:  ADZZ 

RLOUENCV:  1600.0  Hz.  GRAZING  ANGLL: 


PROVINCE:  AEZZ 


AO 


RD-A172  821  THE  PRECISION  OF  BOTTOM  LOSS  PREDICTIONS  FROM  R 

PHVSICRL  SEDIMENT  MODEL (U)  PLRNNIN6  SVSTEHS  INC  NCLERN 
VA  E  J  HOLINELLI  ET  RL.  31  RUB  88  TR-291343 
UNCLASSIFIED  N88814-83-C-86S8  F/B  28/1 


APPENDIX  B 


'•ll5  ys  y  g*rr  'a  '.«■  y  w  'a  wcr  ry  iviv  'A  wt  ww.^v.’^v.-nv  r*  .’ty* 


I 


5 


§ 


> 

> 


v 

t * 


£ 

& 


ri  a) 


in 


uj  Wv 

_j 

i  i 

N  <  ^ 

N  o  -• 

y  Z  2 

CD  n  S 

<*.  “7 

at  ?. 

.  •  a  •-’■ 

U 

CJ  . 

Z  x 

>  o 

O  d 

O'  o 
Q_ 


*T 

r-- 

:c 


-  £ 


t£ 

X 


C 

V  < 

o 

z 


L*  cl 

f—  yt 

i  -  T 


r' 


r  « 


I-  r* 


j-S  ~ 

r-  “ 


V  V  V’ 

y'>f 

$£fe$ 

2£& 


UNCLASSIFIED 


anaamnaimiHgaas 


REPORT  DOCUMENTATION  PAGE 


REPORT  NUMBER 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


4.  TITLE  ('and  Jufctll/.; 

The  Precision  of  Bottom  Loss 
Predictions  from  a.  Physical 
Sediment  Model 


7.  AUTHOAfaJ 

Eugene  J.  Molinelli,  Charles  W.  Holland 
Burlie  A.  Brunson,  and  C.  Pete  Council 


B.  RERFORMINO  ORGANIZATION  NAME  ANO  AOORESS 

Planning  Systems  Inc. 

7900  Westpark  Drive,  Suite  600 
McLean,  Virginia  22102 


II.  CONTROLLING  OFFICE  NAME  ANO  AOORESS 

Defense,  Small  Eusiness  Advanced 

Technology  Program,  Naval  Ocean  Research 
&  Development  Activity,  NSTL,  MS  39529 


4.  PERFORMING  ORG.  REPORT  NUMBER 

TR291343 


*.  CONTRACT  OR  GRANT  numBCRTvj 

N00014-83-C-0650 


10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  *  WORK  UNIT  NUMBERS 


I*.  REPORT  OATE 

31  Auaust  1986 


.  MONITORING  AGENCY  NAME  a  AOORESSfff  different  from  Controlling  Office) 

IS.  SECURITY  CLASS,  (ol  i hit  10 part; 

UNCLASSIFIED 

ISa.  DECLASSIFICATION'  DOWNGRADING 

schedule  n/a 

If.  DISTRIBUTION  STATEMENT  (ol  thle  Report) 


Unlimited 


17.  DISTRIBUTION  STATEMENT  (of  the  ebetrect  entered  In  Block  20 ,  If  different  from  Report) 


If.  KEY  WOROS  (Continue  on  reveree  aide  If  neceeeery  end  Identify  by  block  number) 

Ocean,  Bottom  loss,  continental  terrace,  modeling,  ASW,  seafloor, 
geoacoustic  properties,  sediment  properties 


20.  ABSTRACT  (Continue  on  reeerao  aide  If  neceeeery  end  Identify  by  block  number) 

Ocean  acoustic  signals  and  noise  characteristics  can  be  signifi¬ 
cantly  affected  by  interaction  with  the  seafloor.  One  common 
means  of  quantifying  those  effects  in  terms  of  the  sonar  equa¬ 
tions  is  bottom  loss  (BL) .  BL  can  be  measured  directly  or  com¬ 
puted  from  profiles  of  the  geoacoustic  parameters  density,  com- 
pressional  speed  and  compressional  attenuation.  The  geoacoustic 
properties  can  be  derived  from  physical  properties  of  the 


FOIW 

i  jan  n 


EOlTiON  OF  I  MOV  «S  It  OBSOLETE 
S  N  0)0}-  LF*  01  4-  460 1 


UNCLASSIFIED 

s e curitt  classification  of  this  page  rRn  imara*) 


vl-lvv'.-l'l-l'.  ■  v\'Vv\- 


UNCLASSIFIED 


i-mniu 


sediments,  which  are  more  readily  available,  using  a  physical 
sediment  model  (Biot/Stoll  or  PHYSED)  based  on  the  Biot  theory 
of  acoustic  propagation  in  porous  media. 

Here  we  predict  BL  for  sediments  on  the  continental  terrace  and 
evaluate  those  predictions  in  terms  of  their  precision.  We  use 
readily  available  physical  properties  assembled  into  the  data 
base  PHYPROSE,  with  the  PHYSED  model  to  produce  geoacoustic 
profiles  for  the  calculation  of  BL  by  the  computer  program 
REFLEC.  We  divide  continental  terrace  sediments  into  provinces 
based  on  water  depth  ranges  (shelf  =  0  to  200  m,  slope  =  200  to 
2000  m  and  rise  =  2000  to  4000  m)  and  grain  size  classes  (sands, 
muddy  sands,  sandy  muds,  silts,  etc.).  We  find  that  the  BL  values 
at  a  5°  grazing  angle  for  both  100  Hz  and  1600  Hz  waves,  calcu¬ 
lated  for  each  province,  overlap  to  some  extent,  but  that  useful 
separations  of  BL  values  by  province  do  occur.  Some  provinces 
show  significantly  different  mean  BL  values  (at  the  0.5%  level  of 
significance)  and  some  provinces  show  significantly  reduced  var¬ 
iances  (also  at  the  0.5%  level  of  significance). 

These  results  indicate  a  real  improvement  in  our  capability  to 
predict  BL  if  we  base  those  predictions  on  readily  available 
physical  property  data  and  the  PHYSED  model.  This  result  has 
particular  importance  in  the  shallow  water  environment,  since  I 

sediment  physical  properties  are  often  the  only  readily  obtainable 
information  describing  the  seafloor.  Seldom  are  acoustic  bottom 
loss  or  transmission  loss  data  available.  We  strongly  recommend 
the  application  of  this  physical  sediment  approach  in  shallow 
water  to  further  the  development  of  seafloor  interaction  models 
and  data  bases  suitable  for  operational  use. 
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